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Resumen 



$Que son los agujeros negros? 
$Se expande el universo? $Es concavo o convexo? 
$Quienes somos? $De donde venimos? $A donde vamos? 
frEstamos solos en la galaxia o acompanados? 

- Siniestro Total : $Quienes somos? $De donde venimos? $A 
donde vamos? (1984) - 

of cumulos de galaxias ocupan una posicion relevante en la jerarqma 
de la estructura cosmologica por muchas razones diferentes. Estos, ademas 
de ser una de las estructuras limitadas mas grandes en el universo, contienen 
cientos de galaxias y gas caliente que emite rayos X, por lo que pueden ser 
detectando por su alto corrimiento hacia el rojo. Por lo tanto, estos objetos 
nos permiten estudiar la estructura a gran escala del universo, comprobar 
las teorfas de la formacion de estructura asf como extraer informacion cos- 
mologica dificil de evaluar por otros metodos (ver v.g. Borgani y Guzzo, 
2001; Rosati et al. 2002). 

En el marco del modelo de jerarqui'as para la formacion de estructuras a 
escala cosmica, se piensa que los cumulos de galaxias se formar por acumu- 
lacion de unidades mas pequen as (galaxias, grupos, etc.). Despues de una 
epoca de agregacion (que depende del modelo cosmologico) , le sigue un pro- 
ceso de relajacion violenta que tiende a virializar los cumulos produciendo 
sistemas regulares y quasi-esfericos, con suaves perfiles de densidad tipo 
King. Las recientes observaciones a bajo corrimiento al rojo basadas en 
las orientaciones relativas de las subestructuras dentro de los cumulos, y en 
la relacion entre el estado dinamico y el entorno a gran escala (Plionis y 
Basailakos, 2002) apoyan este escenario jerarquico. 

En la pasada decada, debido al incremento de la resolucion espacial en 
las imagenes de rayos X (ROSAT/PSPC & HRI) y a la disponibilidad de las 
camaras de gran campo, muchos de los cumulos que anteriormente se pensa- 
ban que eran regulares presentan en realidad alguna subestructura. Esto ha 
sido evidente gracias a la aparicion de los satelites Chandra (Weissjopf et al. 
2000) y XMM (Jansen et al. 2001), un hecho que podria tener consecuencias 
importantes para las teorfas de formacion de estructuras. Las propiedades 
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fisicas de los cumulos de galaxias, tales como la fraction de cumulos jovenes, 
las funciones de luminosidad y temperatura, la estructura radial de la ma- 
teria oscura y barionica, asi como las relaciones complejas entre todas estas 
magnitudes, constituyen un desaffo para nuestro conocimiento actual sobre 
como esos objetos crecieron desde las fluctuaciones en la densidad primordial 

Durante los ultimos 20 anos, el paradigma de la materia oscura frfa 
(CDM) ideada por Peebles (1982), Blumenthal et al. (1984) and Davis et 
al. (1985) se ha convertido en el modelo estandar que explica la formation y 
evolution de la estructura cosmica. Este supone que el contenido de materia 
de universo esta dominado por una particula masiva debilmente interactu- 
ante, que todavia no ha sido identificada, y cuyos efectos gravitacionales son 
los responsables de la dinamica de las galaxias y los cumulos. 

Un concepto clave en esta cosmogoma es la formation de los halos de la 
materia oscura. Estos sistemas casi en equilibrio de partfculas de materia 
oscura se forman gracias a un colapso gravitacional no lineal que implica 
periodos de agregacion estacionaria de masa, asf como episodios de fusion 
violenta, en los cuales dos halos de masas comparables colisionan para dar 
lugar a un halo mas grande. 

Se piensa que las galaxias y otros objetos luminosos nacen por el enfri- 
amiento y condensation de bariones dentro de pozos de potential, creados 
por lo halos de materia oscura (White and Rees 1978). En los cmulos de 
galaxias, la mayor parte de los bariones se encuentran en un plasma caliente 
difuso e ionizado que constituye el medio intracumular (ICM). 

La evolution dinamica de las galaxias y del gas ICM plantea un problema 
diffcil para los modelos de formation y evolution de cumulos (v.g. Kaiser 
1986). Aunque la fi'sica asociada a los procesos hidrodinamicos (c.f. Sarazin 
1988) es mas complicada que un colapso gravitacional, los modelos sim- 
ples pueden ser aplicados para relacionar el estado de los bariones con las 
condiciones particulares de la distribution subyacente de la materia oscura.. 
Incluso, se han encontrado algunas correlaciones en las propiedades de los 
rayos X en el medio intracumular, que sugiere que la estructura de densi- 
dad y temperatura de los cumulos puede ser predicha por modelos teoricos 
autoconsistentes. 

El presente trabajo intenta desvelar algunos de los mecanismos que dan 
lugar a las propiedades observadas en los cumulos de galaxias. Nos hemos 
acercado al problema desde diferentes puntos de vista haciendo uso de aprox- 
imaciones anahticas y simulaciones numericas, con el fin de comprender la 
compleja interaction entre las fuerzas gravitacionales a gran escala, que gufa 
la formation de estructura, y la fi'sica barionica, que domina la hidrodinamica 
sobre escalas locales. 

En el capi'tulo uno describimos brevemente el modelo cosmologico asum- 
ido en nuestro estudio. La section 1.1.1. ofrece una somera description 
de los conceptos matematicos basicos necesarios para definir el escenario 
cosmologico. La section 1.1.2 se centra en los valores concretos de los 
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parametros que determinan la geometna del universo, su composition y 
su evolution temporal. La principal evidencia observacional que apoya el 
modelo ACDM tambien sera descrita. 

Hemos hecho una descripcion detallada de los experimented numericos en 
el capftulo 2. Nuestras muestras de los cumulos simulados son el resultado 
de tres codigos independientes basados en muy diferentes escenarios. La 
implementation de los procesos ffsicos en estos algoritmos es brevemente 
resumida, asi'como el procedimiento de analisis que ha sido seguido para 
obtener los resultados que se muestran. Finalmente damos una descripcion 
de las muestras y explicamos los problemas que aparecen en cada uno de 
ellos. 

El capftulo 3 esta dedicado a la distribution de masas que esperamos 
encontrar en los cumulos de galaxias. Especial incapie hacemos en la es- 
tructura radial de los halos de materia oscura, asf como en su origen ffsico. 
Para ser mas preciso: nos centraremos en los perfiles de densidad de nues- 
tras simulaciones de los halos y en su dependencia con el estado dinamico. 
La densidad en el espacio de fases sera tambien investigada, comparando 
los resultados numericos con trabajos anteriores. El origen fisico de la es- 
tructura de la materia oscura en los cumulos sera tratado en el marco del 
formalismo de colapso esferico (Gunn and Gott, 1972; Gunn 1977). Vere- 
mos los perfiles predichos por el modelo del colapso esferico que proporciona 
una descripcion fiel de los experimentos numericos realizados. Intentaremos 
entender las razones por las cuales se producen tal acuerdo, a pesar del he- 
cho de que el crecimiento de los halos de materia oscura siga una serie de 
fusiones no lineales. 

Las propiedades fi'sicas del ICM seran investigadas en el capi'tulo 4. 
En primer lugar, comprobaremos la exactitud de las distintas implementa- 
ciones de las hidrodinamica del gas. Luego testearemos la validez de que 
los cumulos de galaxias estan en un equilibrio hidrodinamicos soportado 
termicamente, asf como en una ecuacion de estado apropiada (politropica o 
isotermica) que permite describir mejor la componente de gas difuso. De- 
spues, intentaremos predecir la densidad radial y los perfiles de temper- 
atura en nuestros simulaciones de cumulos, relacionando estas con las for- 
mas fenomenologicas simples (v.g. Navarro et al. 1997; Moore et al. 1999) 
propuestas para describir los halos de materia oscura. La relation de escala 
entre las propiedades globales de los cumulos, observados en diferentes den- 
sidades, son detalladas en una section diferente. Intentaremos entender el 
origen de la relation masa-temperatura, asi como analizaremos la discrep- 
ancia entre la relation de autosimilaridad en la luminosidad-temperatura y 
la que se observa (v.g. David et al. 1993; Ponman et al. 1996). 

La formation estelar y los efectos de realimentacion van a ser estudiados 
en el capi'tulo 5. Simulamos la historia de la formation estelar cosmologica 
para una serie de escenarios cosmologicos y los comparamos con las observa- 
tions a diferentes frecuencias. Evaluamos los efectos de esta realimentacion 
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en las explosiones de supernovas y la foto-ionizacion de las estrellas jovenes, 
asf como la variancia cosmica inducida cuando un pequenrio volumen es 
considerado. Finalmente, tratamos la cuestion de como el entorno altera la 
actividad de la formacion estelar en las galaxias individuales. La historia 
de la formacion estelar de los cumulos de galaxias es comparada con las de 
camp os de galaxias, e incluimos una breve discusion de los mecanismos re- 
sponsables para de la baja formacion de estrellas halladas en esos ambientes 
tan densos. 

Las principales conclusiones que pueden ser extrafdas de nuestro resul- 
tados son expuestas en la seccion 6. El apendice A contiene un atlas de 
cumulos. Este ensena las caracteristicas individuales de los 15 cumulos que 
reunen la principal caracterfstica analizada en este trabajo. Tambien dibu- 
jamos algunas imagenes similares a las observadas, en las que la diferencia 
entre la cantidad de estructura en gas en materia oscura se muestra con clar- 
idad, asi como el conjunto de figuras en el que las predicciones discutidas en 
el texto principal son aplicadas en una base de cumulo a cumulo. 
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Chapter 1 

Introduction 



Exit light 
Enter night 
Take my hand 

We 're off to never-never land! 

- Metallica : Enter Sandman (1991) - 



clusters occupy a special position in the hierarchy of cosmic 
structure in many respects. Being the largest bound structures in the uni- 
verse, they contain hundreds of galaxies and hot X-ray emitting gas and thus 
can be detected at large redshifts. Therefore, they appear to be ideal tools 
for studying large-scale structure, testing theories of stru cture formation and 
extracting invaluable cq smological information (see e.g. Borgani and Guzzol 



20011: iRosati et al.ll2002h . 



In the framework of the hierarchical model for the formation of cos- 
mic structures, galaxy clusters are supposed to form by accretion of smaller 
units (galaxies, groups etc). After the epoch of aggregation (which depends 
on the cosmological model), violent relaxation processes will tend to viri- 
alise the clusters producing 'regular', quasi-spherical systems, with smooth 
King-like density profiles. Recent low — z obser vations based on the rela- 
tive orientation of substructures within clusters ( West et al. 19951 ) and on 
the relation between their d ynamical state and the large-scale environment 
((Plionis and Basilakosll2~002h do support the hierarchical scenario. 

In the last decade, due to the increased spatial resolution in X-ray imag- 
ing (ROSAT/PSPC & HRI) and to the availability of wide- field cameras, 
many of the previously thought 'regular' clusters ha ye shown to be clumpy 



to some lev el. This is eyen mo re so in the Chandra ( Weisskopf et al. 2000) 



and XMM (j.Tansen et al J l200lh fact that could have important conse- 

quences for structure formation theories. The physical properties of galaxy 
clusters, such as the fraction of dynamically young clusters, the luminosity 
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CHAPTER 1. INTRODUCTION 



and temperature functions, the radial structure of both dark and baryonic 
components, as well as the complex relationships between all these quanti- 
ties, constitute a challenging test for our current understanding of how these 
objects grow from primordial density fluctuations. 

Ove r the la s t 20 years, the Cold Dark Mat ter (CDM) paradig m set out by 
Peebles! (|l982l ). lBlumenthal et ail (|l 9841 ) and bavis etaD (|l98fil ) has become 
the standard model to explain the formation and evolution of cosmic struc- 
tures. It assumes that the matter content of the universe is dominated by an 
as yet unidentified, weakly interacting massive particle, whose gravitational 
effects are responsible for the dynamics of both galaxies and clusters. 

A key concept in this cosmogony is the build-up of dark matter haloes. 
These quasi-equilibrium systems of dark matter particles are formed through 
non-linear gravitational collapse, involving both periods of steady accretion 
of mass, as well as violent merging episodes, in which two haloes of compa- 
rable mass collide to give birth to a larger halo. 

Galaxies and other luminous objects are assumed to form by cooling and 
condensation of baryons within the gravitationa l potential wells created by 
the dark matter haloes ( White and Eeesl 19781 ). In galaxy clusters, most 
baryons can be found in a hot, diffuse, ionised plasma that constitutes the 
intracluster medium (ICM). 

The dynamical evolution of member galaxies and the ICM gas poses a 
diffic ult problem for models of cluster formation and evolution (e.g. iKaiser 



1986). Ne vertheless, alth ough the physics associated to gasdynamical pro- 



cesses (c.f.[^2&) is more complicated than pure gravitational col- 
lapse, simple models can be applied in order to relate the state of the baryons 
to the particular conditions of the underlying dark matter distribution. In- 
deed, several correlations have been observed in the X-ray properties of the 
intracluster medium, hinting that the density and temperature structure of 
clusters may be predicted by self-consistent theoretical models. 

The present work attempts to unveil some of the mechanisms that give 
rise to the observed properties of galaxy clusters. We have approached the 
problem from different points of view, resorting to analytical approximations 
and direct numerical simulations in order to gain some understanding of 
the complex interplay between the large scale gravitational forces driving 
structure formation and the baryonic physics that dominates gas dynamics 
on local scales. 

In this chapter, we briefly describe the cosmological model that will be 
assumed in our study. Section II. 1.11 offers an overview of the basic math- 
ematical concepts required to define a cosmological scenario. Section [1.1.21 
focuses on the actual values of the parameters that determine the geome- 
try of the universe, its composition, and its temporal evolution. The main 
observational evidence supporting the ACDM model will be succinctly de- 
scribed. 

A thorough discussion of the numerical experiments that have been per- 
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formed is given in Chapter [21 Our samples of simulated clusters are the 
outcome of three independent codes, based on very different integration 
schemes. The implementation of physical processes within these algorithms 
is briefly summarised, as well as the analysis procedure that has been fol- 
lowed in order to obtain the results reported here. Finally, a description 
of the samples themselves is given, explaining the problems that are to be 
addressed within each one of them. 

Chapter |3] is devoted to the mass distribution that we expect to find 
in clusters of galaxies. Special emphasis is put in the radial structure of 
the dark matter haloes and its physical origin. To be more precise, we will 
study the density profiles of our simulated haloes and the dependence on 
their dynamical state. The phase-space density will also be investigated, 
comparing the numerical results with previous work. The physical origin of 
the structure of dark matter in clust ers will be addressed within the fr ame- 
work of spherical collapse formalism ( Gunn and Gottlll972l : GunrJl977 ) . We 
will show that the profiles predicted by the spherical infall model provide an 
accurate description of those found in the numerical experiments. We will 
try to understand the reasons for such an agreement, albeit the fact that 
the growth of numerical dark matter haloes proceeds through a non-linear 
series of mergers. 

The physical properties of the ICM will be investigated in Chapter |1] 
First of all, the reliability of current numerical implementations of gasdy- 
namics will be assessed. Then, we will test the validity of the assumption 
that clusters of galaxies are in thermally-supported hydrostatic equilibrium, 
as well as the appropriate equation of state (polytropic or isothermal) that 
better describes the diffuse gas component. Then, we try to predict the ra- 
dial density and temperature profiles of our simulated clusters, relating them 
to the simple phenomenological forms (e.g. iNavarro et al . 1997; iMoore et al 



1999bj) proposed to describe the dark matter haloes. The scaling relation 



between the global properties of the clusters, observed at different overden- 
sities, are dealt with in a separate section. We try to understand the origin 
of the mass-temperature relation, as well as of the reported discrepancy be- 
twee n the self-similar luminosity-temperatu re relation and the observed one 
(e.g. foavid et alJll99.4 IPonman et aflll99fih . 

Star formation and feedback are studied in Chapter 03 The cosmic star 
formation history is computed for a variety of cosmological scenarios and 
compared with observational estimates based on different wavelengths. The 
effects of feedback from supernova explosions and photoionisation by young 
stars are evaluated, as well as the cosmic variance induced when a small 
volume is considered. Finally, we address the question of how does the 
environment alter the star formation activity of individual galaxies. The 
star formation history of clusters of galaxies is compared to that of field 
galaxies, and a brief discussion on the mechanisms responsible for the lower 
star formation found in dense environments is also included. 
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The main conclusions that can be drawn from our results are highlighted 
in Sectional We also quote some of the aspects that in our opinion deserve a 
much deeper investigation than it is possible within the scope of the present 
thesis. 

Appendix lA"l contains a cluster atlas. It shows the individual characteris- 
tics of the 15 clusters that comprise the main sample analysed in this work. 
Mock observational images are plotted, in which the difference between the 
amount of structure in gas and dark matter is clearly seen, as well as a set 
of figures in which the prescriptions discussed in the main text are applied 
in a cluster-to-cluster basis. 

1.1 The cosmological model 

This work (as any other), relies on a series of axioms that it does not intend 
to verify. In our case, maybe the most important of them is the cosmological 
model that is assumed to describe the universe beyond the region of influence 
of our clusters of galaxies. 

The consequences of the chosen cosmological scenario on the formation 
and evolution of galaxy clusters are twofold: On one hand, it sets the power 
spectrum of primordial fluctuations that give rise to cosmological structure. 
On the other, the expansion rate of the universe is largely determined by 
the contribution of its different components to the total energy density. The 
expansion factor, in turn, controls the density of cosmological objects, their 
accretion rate, or the probability that they undergo a merging event. 

1.1.1 Theory 

The basic tenet that governs our current view of the universe is the cosmo- 
logical principle, which is a generalisation of the Copernican idea that the 
Earth does not necessarily have to play a central role in the solar system. 
In general terms, the cosmological principle states that the universe is ho- 
mogeneous and isotropic on large scales, which seems to be corroborated by 
observations of the Large Scale Structure (LSS) or the Cosmic Microwave 
Background (CMB). Furthermore, we will assume that the evolution of the 
universe on such large scales is dominated by gravity, and therefore it can 
be described by the General Theory of Relativity. 

Although it is evident that the universe is anything but homogeneous on 
small scales, the usual approach consists in separating the problem into two 
components: Cosmic structures are considered to grow as a perturbation 
field superimposed to an otherwise homogeneous and isotropic space-time 
fabric. The evolution of the universe as a whole is reduced, in this simple ap- 
proximation, to the dynamics of cosmological expansion, which constitutes 
one of the most fundamental applications of General Relativity. 
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Friedmann equations 

It can be shown (e.g. that the assumptions of homogeneity 

and isotropy are equivalent to the requirement that the metric tensor of the 
universe takes the form of the Robertson- Walker metric 



ds 2 = dt 2 - a{tf 



— ? + r 2 {d6 2 + sin0#) 

1 — kr z 



[1.1) 



= R^u - t; r 9^ = 87rGT^ + Ag^ (1.2) 



where the dimensionless function a(t) is the cosmic expansion factor. 

According to General Relativity, the geometry of the space-time back- 
ground is related to its matter-energy content by virtue of the Einstein field 
equations, 

1 

2' 

The stress-energy tensor of the universe is taken to be, also from the 
cosmological principle, that of a perfect fluid (such that a comoving observer 
would observe an isotropic universe), 

T^v = P9fiu + (p + p)u^u u (1.3) 

Inserting the perfect fluid tensor and the Robertson- Walker metric into 
the Einstein field equations above, we can apply the stress-energy conser- 
vation law {T^ u = 0) to get the Friedmann equations, which govern the 
dynamics of the cosmic scale factor a(t): 

d\ 2 8vrG A k d 4vrG . A 
a J 6 6 a z a 6 6 

The solutions of these equations are known as Lemaitre-Friedmann- 
Robertson- Walker (LFRW) universes, and their characteristic properties 
vary greatly as a function of the several densities that can appear in (|1.4|) 
and their corresponding equations of state. 



Density parameters 

Before the 1990's, little could be said about the matter-energy composition 
of the universe. The only points in which a general agreement existed were 
the validity of the cosmological principle, the finiteness of the age of the 
universe, and its expansion from a gas cloud of extremely high density and 
temperature. 

This 'hot Big-Bang' scenario had been induced and confirmed by several 
independent obs ervational evi dences, namely the measurement of redshift in 
nearby galaxies the relative abundan ces of light elem ents (ex- 

plained through a primordial nucle osynthesis phase, Gamow 1946), and the 

very 



detection of the CMB radiation by Penzias and Wilsonl (|l965l ). With 
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Density par. State eq. 



Matter Relativistic f2 7 + Ohdm p = 1/3/9 

Non-relativistic + ^cdm p = 

Space-time Vacuum P = —p 

Curvature 1\ P = — 1/3/9 



Table 1.1: Possible energy densities in the universe 



few ingredients, this simple model was able to explain satisfactorily enough 
the main features of the observable universe from its very first seconds of 
existence. 

Nevertheless, LFRW universes have a set of free parameters that only in 
the last decade has it been possible to bound. On one hand, the pressure 
and density terms that appear in (|l,4j) corresponding to the different types 
of energy density that may exist in the universe, and, on the other hand, 
the actual value of the expansion rate at the present time, usually known as 
Hubble constant and expressed as 

Hit) = ^r\ ; H = H(t ) = 100 h (km s^Mpc" 1 ) (1.5) 
a(t) 

It is useful, in order to classify the different terms contributing to the 
stress-energy tensor, to work with the density parameters, defined as frac- 
tions of the critical density required for the universe to have a spatially-flat 
geometry (k = 0): 

In analogy with ordinary matter, we define an effective density due to 
the curvature of space-time and to the cosmological constant term (which 
can be interpreted as a uniform vacuum energy density p\ = £l\p c ). 

- ^ : «a(*) = (1.7) 

With this definitions, the first Friedmann equation transforms into the 
following closure relation: 

£>,(*) = 1 (1.8) 

i 

Table 11.11 summarises a classification of the possible components that 
may contribute to the total energy density of the universe, according to the 
equation of state relating their density an pressure. These equations also 
define the relationship between the different densities and the cosmic scale 
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Figure 1.1: Cosmic triangle 



factor a(t). To obtain the time-dependence f2j(a) we just have to insert the 
corresponding equation of state in (|1.4I) : 

Pi (t) = u iPi (t) Pi (t) = pi(to) a- 3( - 1+ ^ (1.9) 

As a result, the current values of the density parameters completely 
determine the solutions of the Friedmann equations. From now on, if there 
is no explicit time dependence of any density, we will take 0^ to denote the 
present-time value (unless otherwise noted). The first Friedmann equation 
may then be written as 

H \ 2 

— = Q r a" 4 + n m a~ 3 + n A + n k a~ 2 (1.10) 
Ho J 

where it becomes evident that the properties and evolution of a homogeneous 
and isotropic universe depend only on the four density parameters (subject 

to the closure relation ll.8j) plus the Hubble constant. 

This is very aesthetically summarised in Figure rmrtBahcall et al.lll999h . 
where the parameter space is represented by the position of the model uni- 
verse in the 'cosmic triangle'. Each point within the triangle satisfies the 
closure relation l|1.8|) . The horizontal line (marked 'flat') corresponds to 
k = 0, separating an open universe from a closed one. The red line, nearly 
along the A = line, separates a universe that will expand forever (approx- 
imately A > 0) from other that will eventually re-collapse (approximately 
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A < 0). And the yellow, nearly vertical line separates a universe with an 
expansion rate that is currently decelerating from one that is accelerating. 

The locations of three key models are highlighted in the figure: The 
Standard Cold Dark Matter (SCDM) model, dominated by non-relativistic 
non-baryonic matter, (O m , ^a) = (1, 0); an open universe with no cosmolog- 
ical constant, (Sl m , I^a) = (0.3, 0); and a flat universe with a A-term, defined 
by {n m ,n A ) = (0.3,0.7). 



1.1.2 Observations 

From a theoretical point of view, there are no preferred values of the density 
parameters. Therefore, constraints on their actual values in the real uni- 
verse must come from observational estimates. The last ten to fifteen years 
have witnessed a tremendous improvement in observational instrumenta- 
tion. These novel technologies have opened a plethora of new possibilities 
to measure the contribution of the different types of energy densities to the 
total cosmological budget. 

While in the 1960s Sandage's goal was to 'search for two numbers', it 
is now appreciated that even the basic cosmological model requires about a 
dozen or so parameters, which are linked in non-trivial ways. Combinations 
of several independent experiments must be used to break these degenera- 
cies. In this section, we review some recent results aimed to find the values 
of the cosmic density parameters that best fit the whole set of available ob- 
servations. As we shall see below, all observational evidence 1 seem to favour 
the so-called ACDM cosmological scenario. 



Baryon density 



The standard theory of Big-Bang Nucleosynthesis (BBN) predicts the abun- 
dances of the light element nuclei H, D, 3 He, 4 He, and 7 Li a s a function of the 
cosm ological baryon-to-photon ratio, r] = n^jn^ (see e.g. iKolb and Turner 
1990). A measurement of the ratio of any two primordial abundances gives 
ij, and hence the baryon density, while a second ratio tests the theory. 

However, it is extremely difficult to measure primordial abunda nces be- 
cause in most places gas ejected from stars has enriched the medium. Adamsl 
(jl97fih suggested that it might be possible to measure the primordial D/H 
ratio in absorptio n -line systems toward QS Os . Several estimation s (e.g. 

have 



Tvtler et al 



- -19961 : buries and Tvtleil Il998al lO'Meara etaA 12001. 
been made thereafter. The baryon density reported bv lO'Meara et a*l*l i2001) 
amounts to 

(1.11) 



tt h h 2 = 0.0205 ± 0.0018 



1 At the moment of writing this thesis. 
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Hubble constant 

In standard Big Bang cosmology, the universe expands uniformly; and lo- 
cally, according to the Hubble law, v = H^d, where v is the recession velocity 
of a galaxy at a distance d, and Hq is the Hubble constant, the expansion 
rate at t he cu rrent epoch. More than seven decades have now passed since 
Hubblel (|l929l ) initially published the correlation between the distances to 



galaxies and their recession velocities, thereby providing evidence for the 
expansion of the universe. But pinning down an accurate value for the Hub- 
ble constant has proved extremely challenging. There are many reasons for 
this difficulty, but primary among them is the basic difficulty of establishing 
accurate distances over cosmologically significant scales. 

Measuring an accurate value of Hq was one of the motivating reasons for 
building the NASA/ESA Hubble Space Telescope (HST). The overall goal of 
the Hq Key Project was to measure the Hubble constant based on a Cepheid 
calibration of a number of independent secondary distance determination 
methods. To extend the distance scale beyond the range of the Cepheids, 
a number of methods that provide relative distances were chosen, including 
the Type la supernovae, the Tully-Fisher relation, the fundamental plane for 
elliptical galaxies, surf ace-brightness fluctua tions, and Type II supernovae. 

A recent paper bv lMonld et"aTI <|2000al lbh combines the results of the 



HST with these secondary methods base d on different standar d candles. 
The final results of the Hq Key Project (|Freedman et al. 20011 ) yield the 
following value: 

H = 72 ± 8 km s" 1 Mpc" 1 (1.12) 



Decelration parameter 

Just as large distance measurements on Earth show us the curvature on 
Earth's surface, so do large distance measurements in cosmology show us 
the geometry of the universe. Since the geometry of space-time is given by 
the values of the cosmic density parameters, measurements of the luminosity 
distance at high redshifts can be very help to constrain those values. 

In these sense, observatio ns of Type la supernov ae at z > 3 by the 



in tnese sense, observatio ns ot Type la supernovae at z > 6 by trie 
Supernova Cosmology Project (Perlmutter et al.lll999h and the High-Z Su- 



pernova Search Team IjRiess et al.lll998l : ISchmidt et al.lll998l ) have provided 
an extremely useful benchmark to test the different cosmological scenarios. 
These measurements allowed for the first time an accurate determination 
of the deceleration parameter qQ = —ad/ a 2 = Q m /2 — S1a- The most strik- 
ing conclusion drawn from these experiments was that the expansion of the 
universe is nowadays accelerating (qQ < 0), which hints the presence of a 
non-zero positive cosmological constant, 



(1.13) 
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Curvature 



One of the most promising techniques to determine the geometry of the 
universe is the study of the amplitude fluctuations in the cosmic microwave 
background. Early detection of a p eak in the region of the so-called first 
acoustic peak ([Netterfield et al 1ll99flL as well as the availability o f fast codes 
to compute theoretical amplitudes ( Sefiak and Za ldarriaga 1996) , have pro- 
vided a first constraint on the value of (|Lineweaver et al.ll997l : Hancock et alJ 
Il998h . 

The spectacular results of Boomerang and Maxima balloon experiments 
have firmly established that geometry of the universe is very close to flat 
dde Bernardis et a1.ll2000l:lHauanv et aljEoOfl: iBalbi et aljEoOrl: ll.anPO el al. 
200 If ) . Recently, iBenoit et alJ (|2002l ) combined t he anisotropics measured by 



the Archeops experimen t with d ata from C OBE (jTegrnarlJl99fih . B oomerang 
dNetterfield et al.ll2002l). Dasi (jrlalverson et al.ll2002h . Maxima dLee et al. 



200lh . VSA (jScott et al1l2002T ) and CB1 (jPearson et al .1 Eool ) . The com 



bination of all these CMB experiments provides an estimation of the total 
energy density of the universe 



1 04+ 012 



(1.14) 



as well as a measure of the baryon density completely independent (but 
consistent) with the nucleosynthesis estimate: 



n h h 2 = 0.022" 



-0.003 
-0.004 



(1.15) 



Adding other constraints, such as those imposed by the observed large 
scale st ructure (see e.g. the an alysis based on the 2dF Galaxy Redshit Sur- 
vey by lEfstathiou et all 20021) or the o bserved baryon fraction in galaxy 
clusters (e.g. ISadat and Blanchardl 200 ll ) , the values of the density parame- 
ters are still compatible with those listed above. The so-called 'concordance' 
model emerges as an apparently solid observational measurement, confirmed 
by several independent methods. 

Given the amount of observational evidence that currently supports the 
ACDM cosmology, we will base most of our study on the assumption that 
the space-time background is well described by a homogeneous and isotropic 
LFRW universe with ft m = 0.3, £l A = 0.7, tt h = 0.04 and h = 0.7. In those 
cases where some other cosmology is used, it will be explicitly noted. 



Chapter 2 

Numerical experiments 



One more, one more, one more, one more... 
Addict with an Apple Mac, megabyte maniac 
Moron with a mouse mat, a junky 
Distressed, you've guessed, obviously obsessesed 
I, NEED MORE MEMORY! 

- Toy Dolls : One More Megabyte (1997) - 



rllllietical simulations of three-dimensional self-gravitating fluids have 
become an indispensable tool in cosmology. They are now routinely used to 
study the non-linear gravitational clustering of dark matter, the formation 
of galaxy clusters, the interactions of isolated galaxies, and the evolution of 
the intergalactic gas. Without numerical techniques the immense progress 
made in these fields would have been nearly impossible, since analytic cal- 
culations are often restricted to idealised problems of high symmetry, or to 
approximate treatments of inherently non-linear problems. 

The advances in numerical simulations have become possible both by 
the rapid growth of computer performance and by the implementation of 
ever more sophisticated numerical algorithms. The development of powerful 
simulation codes, with high level of parallelism, still remains a primary task 
in order to take full adva ntage of new computer technologies. 

Early simulations fe.g. lWhitell97dlFalil97^ : lAarseth et al.ll979l ) largely 
employed the direct summation method for the gravitational N-body prob- 
lem, which remains useful in collisional stellar dynamical systems, but it is 
inefficient for cosmological simulations with larg e iV due to the rapid in- 
crease of its computational cost with ./V (see e.g. for a review 
on numerical methods for cosmological simulations). 

In order to overcome the problems inherent to the direct summation 
N-body methods, a large number of groups have developed new techniques 
for collisionless dynamics that compute the large-scale gravitational field 
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in a (regular or irregular) grid. In the simplest implementation of grid- 
based N-body methods, the P article-Mesh (PM) N-body technique (see e.g. 
Hocknev and Eastwood! Il98lh , the Poisson equation is solved using Fast 



Fourier Transform (FFT) in a regular grid of the density field, that is con- 
structed by interpolation from the particle positions to the grid nodes. New- 
tonian forces are computed by differentiating the potential on the grid and 
interpolating back to the particle positions. The time for force computa- 
tion scales as 0{N log N) in these methods, which makes them able to treat 
a large amount of particles. In fact, this is the N-body method that can 
handle the largest number of particles, and it has been used widely to in- 
vestigate the non -linear clustering of different cosmological mode ls since the 
early 1980's fe.g. lKlvpin and Shandarin1ll98.4 Iwhite et alJll98.^ . 



The lack of resolution belo w the 
the PM model. iHocknev et al. 



$rid-size is one of the major problems of 
1) proposed an hybrid scheme to increase 
the numerical resolution of the PM method, consisting in the decomposition 
of the force acting on each particle into a long-range force, computed by the 
particle-mesh interaction, and a short-range force due to nearby particles. 
Modern versions of these codes use finer grids on highly clustered regions to 
recursively apply the P 3 M algorithm. Therefore, the number of neighbours 
per particle is considerably reduced, which makes the code much faster with 
respect to the original implemen tation. Th e publ icly available Adaptive 
P 3 M (AP 3 M) code developed by Couchmanl ( 1991 ) is a good example of 
this technique and has helped the development of other codes based on this 
integration scheme. 

Another approach to increase resolution of PM methods is to use non- 
uniform grids and algorithms to adapt the computational mesh to the struc- 
tures formed by gravitational clustering. The Poisson equation can be solved 
on a hierarchically refined mesh by means of finite-difference relaxation 
m ethods, an approach ta ken in the Adaptive Refinement Tree ( ART) code 



bv lKravtsov et alJ (1199711 . the MLAMP code (jKnebe et al 
recently bv lTevssierl fcool ) in the RAMSES AMR code. 



2001), and more 



An alt ernative to these schemes are the so-called tree algorithms, pio- 
neered bv lColella and Woodward! (|l985h . Tree algorithms arrange particles 
in a hierarchy of groups, and compute the gravitational field at a given 
point by summing over multipole expansions of these groups. In this way 
the computational cost of a complete force evaluation can be reduced to a 
0(N log N) scaling. The grouping itself can be achieved in various ways, for 
example with Eulerian sub divisions of space ( Barnes and Hutlll98(jh. or w ith 
nearest-neighbour pairings flPress et al.lll98 d: ljernigan and Porterl ll989\ A 
tech nique related to o rdinary tree algorithms is the fast multipole-method 
(e.g. iGreengardl 1988), where multipole expansions are carried out for the 
gravitational field in a region of space. 

While mesh-based codes are generally much faster for nearly homoge- 
neous particle distributions, tree codes can adapt flexibly to any clustering 
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state without significant losses in speed. This Lagrangian nature is a great 
advantage if a large dynamic range in density needs to be covered. Here 
tree codes can outperform mesh based algorithms. In addition, tree codes 
are basically free from any geometrical restrictions, and they can be eas- 
ily combined with integration schemes that advance particles on individual 
timesteps. 

An independent way to speed up force computations for large number 
of particles, is to use efficient parallel algorithms to distribute the work load 
among many processors. A big effort has been done rece ntly in develor 



among many processors. A Dig enort nas been done rece ntly m develop- 
ing p arallel N-body cod es, usin g either direct summatio n rtMakind f2002). 



PM (iRicker et al.ll200(tl. P 3 M (MacfarlanH etabl ll 998F) or tree methods 
(|Dubinskilll996l : iLia and Carraro T2OO1I : ISpringel et al, 2001bh . The advan- 



tage of parallel N-body codes is the portability. They can run from cheap 
Beowulf PC clusters to very expensive supercomputers. 

In recent years, collisionless dynamics has also been coupled to gas dy- 
namics, allowing a more direct link to observable quantities. Traditionally, 
hydro dynamical simulations have usually employed some kind of (Eulerian 
or Lagrangian) mesh to represent the dynamical quantities of the fluid. Mod- 
ern grid-based computational fluid dynamical codes are based on Godunov 
methods: i.e. gasdynamical equations are integrated in each volume ele- 
ment of the computational mesh. The evolution of volume- aver aged fluid 
quantities can be computed from fluxes through the cell boundaries, which 
are estimated from approximate Riemann solvers. 

A particular strength of these codes is their ability to accurately resolve 
shocks without artificial viscosity terms that might introduce numerical dis- 
sipation. A regular mesh also imposes restrictions on the geometry of the 
problem, as well as onto the dynamic range of spatial scales that can be 
simulated. New adaptive mesh refin ement codes (fNorman and Brvanl ll999: 
Kravtsov et al.ll2002l : lTevssierll2flf)2l ) have been developed to provide a solu- 
tion to this problem. 

In cosmological applications, it is often sufficient to d escribe the g as by 
smoothed particle hydrodynam ics (SPH) , as invented by iLucvl ( 1977T ) and 
lOing old ht,H Mo^harJ » The particle-based SPH iseldremely flexi- 
ble in its ability to adapt to any given geometry. Moreover, its Lagrangian 
nature allows a locally changing resolution that 'automatically' follows the 
local gas density. This convenient feature helps to save computing time by 
focusing the computational effort on those regions that have the largest gas 
concentrations. Furthermore, SPH ties naturally into the N-body approach 
for self-gravity, and can be easily implemented in three dimensions. 

These advantages have led a number of authors to develop SPH codes 
for applications in cosmology in combination with different flavours of N- 
body codes, to the point that nowadays these hydrodynamic codes are by 
far the most numerous in cosm ological simulations. A comparison study can 



be found in iFrenk et al. I (jl999h . where the major differences between SPH 
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and grid-based gasdynamical codes were found in the central part of the 
radial entropy and temperature profiles of the galaxy cluster used for the 
comparison. The AMR gasdynamical codes tend to produce an isentropic 
gas profile in the inner regions, while SPH codes predict an isothermal profile 
with a decreasing entropy towards the centre of cluster-size haloes. We will 
come back to this problem in Section 14.11 showing that results from SPH 
and grid-based codes can be reconciled wh en an appropriate treatment of 
entropy conserva tion in SPH codes is used ((Springel and Hernauisl]l2002al : 
ISerna et al .Il2002h . 

The present chapter is devoted to the description of the numerical ex- 
periments reported in this thesis. First, the different codes that have been 
used are described in some detail, focusing on the specific merits of each one 
depending on the sort of problem that we wanted to tackle. The general 
procedure followed to analyse the results is discussed on a separate section, 
and finally a description of the actual numerical experiments is given. 



2.1 Description of the codes 

In order to study the internal structure of dark matter haloes, a high- 
performance N-body algorithm is required to accurately resolve their very 
central parts. For this purpose we used two independent state-of-the-art 
N-body codes; the parallel Tree-SPH code Gadget and the pure N-body 
implementation of the AMR code ART. Results based on these two inte- 
gration schemes have been compared in order to assess the reliability of our 
conclusions (see Section EOI) . 

Inclusion of gas dynamics in N-body simulations is even more demanding 
from the computational point of view, enforcing the use of efficient paral- 
lel codes for this purpose. The dynamical evolution of the gas component 
within dark matter haloes has been investigated by means of a series of 
simulations accomplished with Gadget. Because of the problems with en- 
tropy conservation in SPH codes, a new implementation of Gadget has 
been used in which the entropy (rather than energy) conservation equation 
is integrated ( Springel and Hernquistl 2002al ) . As a test of the consistency 



of our results, we have compared in Section 14.11 one of the galaxy clusters 
simulated by means of the SPH technique with an independent numerical 
experiment by Nagai and Kravtsovl rt2002h usi ng the new hydrodynamical 



version of the APT code (|Kravtsov et al.ll2002l) . 

Finally, to address the problem of star formation in different environ- 
ments, it does not suffice to take into account only gravity and gasdynam- 
ics, but also all the complex physics related to star-gas interactions. In 
this regard, we have used a Particle-Mesh code for N-body coupled with an 
Eulerian hydrodynamical code based on the higher order Godunov Piece- 
wise Parabolic Method (PPM). The code is supplemented with a model for 
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cooling, star formation and feedback from supernovae. From now on, we 
will refer to this c ode as YK 3 , which stands for the initials of the authors 
( Yenes et a,l1ll997t V 

All the codes employed in this work have been parallelised to run effi- 
ciently either in distributed (Gadget) or Shared Memory Processor systems 
(ART and YK 3 ). Numerical experiments have been performed in a variety 
of supercomputers and centres: 



• SGI Origin 2000 at the CEPBA (Barcelona, Spain) 

• SGI Origin 2000 at the NCSA (Chicago, USA) 

• SGI Origin 3800 at the CIEMAT (Madrid, Spain) 

• Hitachi SVR at the AIP (Potsdam, Germany) 

• Hitachi SVR at the LRZ (Munich, Germany) 



2.1.1 YK 3 

This code implements the most complete treatment of the physical processes, 
at the expense of a much lower spatial resolution than that achieved by ART 
or Gadget. 

The matter in the simulated Universe consists of four phases. 1) The 
dark matter (labelled by a subscript "dm" ) in th form of weakly interacting 
collisionless particles is the main contribution to the mean density of the 
universe (Odm = 1 — ^b)- The baryonic component is described as a medium 
consisting of three interacting phases: 2) hot gas (labelled by subscript h, 
T h > 2xl0 4 K), 3) gas in the form of cold dense clouds (subscript c, internal 
temperature T c = 10 K) resulting from cooling of the hot gas, and 4) 
"stars" (subscript *), formed inside cold clouds and treated as collisionless 
particles. Thus, the total density (r) is the sum of four components: 

P = Pdm + Ph + Pc + P*- (2.1) 

This picture is consistent with work on models of galaxy formation and 
evolution (see for example Nulsen 1986; Thomas, 1988a,b; Hensler and Burk- 
ert, 1990; Daines et al 1994; Nulsen Fabian, 1995) and in this context 
appears to be superior to a treatment of the gas component as a one-phase 
medium. 



Gravity and Gas Dynamics 

Gravity is described by the standard Particle-Mesh algorithm, i.e. density is 
computed at the grid points by Cloud in Cell interpolation from the particle 
positions. Poisson's equation is the solved on the grid by FFT and forces are 
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obtained by differentiating the potential in the grid. Forces acting on each 
particle are found by interpolation from the grid to the particle positions 
using the same kernel (CIC). Finally, Newton equations are integrated for 
each particle using a Leap-frog scheme, to advance them. 

Gas dynamics is described by the high order Godunov PPM algorithm 
( Colella and Woodward! 1984 ). As we said at the beginning of this section, 
the main advantage of this Eulerian code is that it does not involve an 
artificial viscosity, but still treats shocks very accurately. This method is 
very fast and highly parallelisable. This algorithm is applied to the gas- 
dynamical equations in one dimension at a time. Multiple dimensions are 
treated by directional timestep splitting, whereas local processes (heating, 
cooling, star formation, etc.), which involve various components of collisional 
matter, as well as self-gravitation and gravitational interaction with dark 
matter, are treated using process timestep splitting (Oran & Boris 1986). 



Additional physics 

A proper treatment of cooling and star formation requires consideration of 
processes occurring on scales well below the numerical resolution of the code. 
Despite certain simplifications, the approach taken in YK 3 relies heavily on 
the picture of the interstellar medium by McKee & Ostriker (1977), who 
provided a detailed description of the physics in dense regions. This theory 
assumes the existence of cold clouds forming due to thermal instability in 
approximate pressure equilibrium with the surrounding hot gas. 

Despite the existence of multiple phases with different temperatures, the 
cooling rate can still be computed from the mean parameters of the hot gas. 
Let us denote the local cooling rate of the plasma due to radiative processes 
by dE/dt = —A r (p,T), where p and T here represent the true local values 
of the gas density and temperature. According to the theory of McKee and 
Ostriker, the rate of energy loss expressed in terms of average gas density 
Ph and temperature T^ will be higher than the nominal rate — A r (p/ l , T^) by 
a cooling enhancement factor C, which is taken to be C = 10 to account for 
all the effects resulting from unresolved density inhomogeneities. 

Cooling depends also on the chemical composition of the gas. In order 
to incorporate the effects of metallicity into the code, solar abundance is 
assumed in regions where either thermal instability is present or previous 
star formation has occurred. Otherwise, the region is considered not to be 
enriched by metals, and the cooling rate for a gas of primordial composition 
(Fall & Rees 1985) is used. 

The ionisation of the intergalactic and interstellar medium by UV pho- 
tons emitted by quasars, AGNs, and by nonlinear structures in the process 
of formation has two important consequences: First, gas outside dense re- 
gions (/>gas < 2(/9 gas ), where (p ga s) is the mean cosmological gas density) 
is ionised and has a temperature of T = 10 38 K (Giroux & Shapiro, 1996; 
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Petitjean, Miicket, & Kates, 1996; Miicket et al., 1996). Second, in low to 
moderate-density regions, heating due to the ionising flux suppresses the 
thermal instability (Miicket & Kates, 1996) and thus the formation of cold 
clouds and ultimately stars. Hence, cold clouds are allowed to form only if 

/Ogas/Pcr > T> • ^bar, with V ~ 50 - 100. 

In addition to the above processes, energy loss due to Compton cooling 
is also included, given by Arjomp = 7 x 10~ 36 n#T e a -4 , where nu is the 
number density of hydrogen ions, T e is the electron temperature, and a is 
the expansion parameter. 

In this model, there are two mechanisms for gas to leave the hot phase 
and enter the cold phase. First, if the temperature of the hot gas drops 
below a threshold temperature Ty mi = 2 x 10 4 K, and p gas /p C r > T> ■ O^ar 
all the hot gas is transferred immediately to the cold phase, thus making it 
available for star formation. This process gives the following terms in the 
continuity equations for the gas: 



( dph \ - Ph t €h pgas n T) o T ^ T 
—JT — ~1 > *cool — q 7^7; > ^'"bar J- < ilim, 

V dt ) hot^coid *cooi de h /dt p CT 

(2.2) 

where i coo i is the cooling time, and eh is the thermal energy per unit mass. 
The transfer may be treated as immediate if the cooling time is very short 
compared to a timestep, which typically happens at T pt 2 x 10 4 K. 

A second way of transferring gas to the cold phase is by sufficiently 
rapid thermal instability. The temperature range for creation of cold clouds 
depends in reality on the ionising flux and the local density (Miicket & 
Kates, 1996), but the simplified restriction T < T; nst = 2 x 10 5 is used in 
the code. To estimate the rate of growth of mass in cold clouds, the energy 
emitted by the hot gas is actually lost. If and e c are the internal thermal 
energies per unit mass of hot and cold gas, then the change of energy of the 
system due to radiative cooling is 

MPheh + Pc£c - =-CA r ( Ph ,T h ), (2.3) 



V dt 



cooling 



where C is the enhanced cooling factor due to unresolved dumpiness. We 
take e c = constant, which means that the cold gas cannot cool below T\ im . 
For an ideal gas, these assumptions imply extra terms in the continuity 
equations: 

dfh \ = _(^Pl\ CA r (p h ,T h ) 

) therm.inst V dt J therm. inst 7 e h ~ £c 

where 7 is the ratio of specific heats. 

Because of the frequent exchange of mass between the hot and cold gas 
phases, it is reasonable to consider the hot gas and cold clouds as one fluid 
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with rather complicated chemical reactions going on within it. Thus, only 
the motion of the hot component is computed in the code, integrating the 
change of the total density of the gas (/9 ga s = Pc + Ph)- 

Both phases are assumed to be in pressure equilibrium, P gas = P c = 
Ph = (l — l)wh) where « n = ph e h is the internal energy per unit volume. 
The velocity associated with the fluid involves an average over a cell size, but 
in a multiphase medium with supernovae, the hot gas is actually "windy" , 
and the cold clouds have a significant velocity dispersion (McKee & Ostriker, 
1977; Cowie, et al., 1981; Hensler and Burkert, 1990.). Therefore, the code 
includes an effective pressure Pdisp °c p c Th associated with the cold gas 
component. 

Star formation takes place in the cold clouds, leading to a decrease in 
the density given by 



where t* ~ 10 8 yr is a fixed characteristic star formation time. The actual 
time-scale for star formation can exceed t* and will depend strongly on the 
rate of the conversion of hot gas to the cold phase. The lifetime of massive 
stars is about 10 7 yr or even shorter; thus most of stars with M* heavier 
than (10 — 20)Mq will explode as supernovae during one timestep. Stars 
in the mass interval from (5-7)Mq to IOMq will explode on a longer time- 
scale, but they produce less energy, and therefore in view of the uncertainties 
in supernova energy, their energy input is neglected in the model. Due to 
the explosion of massive stars, the stellar growth rate is decreased by the 
fraction (3 of stars that explode as supernovae: 



The precise value of (5 is very sensitive to to the form of the initial 
mass function (IMF), particularly its low- mass limit. The principal source 
of uncertainty is the dependence of the IMF on the abundance of heavy 
elements, which can be of critical importance at very early epochs of galaxy 
formation. In view of these uncertainties, (3 is taken to be constant in time. 
For a Salpeter IMF, the fraction of mass in stars with mass larger than 
1OM is = 0.12. 

Evaporation of cold clouds is an important effect of supernovae on the 
interstellar medium (McKee & Ostriker 1977, Lada 1985). The total mass 
of cold gas heated and transferred back to the hot gas phase is assumed 
to be a factor A higher than that of the supernova itself. The supernova 
feedback parameter A could depend on, among other things, the energy of 
the supernovae, the cloud spectrum, and the ambient density. McKee and 
Ostriker (1977) give an estimate of the evaporated mass which scales with 




star— formation 



Pc 



(2.5) 



dp* = (1 - (3)p c 
alt t* 



(2.6) 
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the energy of the supernova E and the gas density as E®^ 5 •n h 4/ ' 5 . However, 
at present we do not include the dependence on the gas density, because 
here refers to the local density, which is relevant for the propagation of 
the supernova shock, whereas our is the mean density of the hot gas 
averaged over quite large a volume. 

Assuming that the energy input due to supernovae is proportional to 
their total mass, we obtain for the evaporation rate 

dph\ _ (dp c \ _A(ip c 



/evap V dt / evap t* 

Correspondingly, the net energy supplied to the hot gas phase by super- 
novae is: 

where e SN = (^sn)/(M sn ), (#sn) = 10 51 erg, and (M SN ) = 22M (Salpeter 
IMF). Admissible values of the supernova feedback parameter A are con- 
strained by the condition that the energy of a supernova explosion must 
be larger than the energy required for the evaporation of cold clouds. This 
restriction leads to 6sn > A ■ (eh — e c ). For the above supernova energy and 
mass and for hot gas temperature Th ~ 10 6 K, this condition implies the 
restriction 1 < A < 250. 

Therefore, for large values of A, cloud evaporation dominates over ther- 
mal re-heating, which translates into less pressure gradients and more mass 
transfer between cold and hot phases. 



2.1.2 ART 

The spatial resolution of particle-mesh codes like YK 3 can be substantially 
improved by means of adaptive grids that provide more resolution in the 
high density regions where galaxies form. The Adaptive Mesh Refinement 
(AMR) method increases the dynamical range with respect to YK 3 by about 
two orders of magnitude without loss of mass resolution or computational 
speed. 



Gravity 

The A daptive Refinement Tree (ART) N-body code developed bv lKravtsov et al 



(1997) is based on this technique. The computational volume is covered by 
a cubic regular grid that defines the minimum resolution, where the Poisson 
equation is solved with a traditional Fast Fourier Transform (FFT) tech- 
nique using periodic boundary conditions, as in the PM algorithm. The 
difference is that this grid is refined (i.e. its cells are divided) wherever the 
density exceeds a predefined threshold. 
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Any cell can be subject to further refinements; the local refinement pro- 
cess continues recursively until the density criterion is satisfied. Once con- 
structed, the mesh is adjusted at each timestep to the evolving particle 
distribution. Finer meshes are built as collections of cubic, non-overlapping 
cells of various sizes organised in octal threaded trees (i.e. each node may be 
split in eight sub-cells). 

This is the same construct used in standard Tree codes, but there is an 
important difference in tha t cells are conne cted to their 6 neighbours in a 
Fully Threaded Tree (FTT, iKhokhlovl f998) on all levels of refinement. In 



addition, cells that belong to different trees are connected to each other 
across tree boundaries. We can consider all cells as belonging to a single 
threaded tree with a root being the entire computational domain and the 
base grid being one of the tree levels. 

The refinement procedure can be used either to construct the mesh hier- 
archy from scratch or to modify the existing meshes. However, in the course 
of a simulation the structure is neither constructed nor destroyed. Instead, 
the existing meshes are modified every computational cycle to account for 
the changes in particle distribution. Therefore, the computational mesh is 
adapted automatically to the structures developed as clustering proceeds. 

Poisson equation must be solved on every refinement level. Therefore, 
the value of the density must be computed for every cell regardless of 
whether or not it is a leaf. On each level, starting from the finest one 
and up to the zeroth level, the density is assigned using the standard CIC 
interpolation from the particle positions. For the zeroth level of the mesh 
hierarchy (regular grid of fixed resolution), Poisson equation is solved by 
the stan dard FFT method. In the refinement sub-mesh es, the relaxation 
method (|Hocknev and Kastwoodl Fl98l1: IPress et alJ lTflSfi) is used. Masses 



and timesteps are independent for each particle. 
Gas dynamics 

The ART code achieves high spatial resolution by adaptively refining regions 
of interest using an automated refinement algorithm. Due to the use of an 
Eulerian mesh hierarchy in the FTT algorithm, the inclusion of Eulerian 
gasdynamics is a natural extension to the pure N-body part. 

As in YK 3 , the equations of g asdynamics and part i cle m o tion are inte- 
grated in 'supercomoving' variables lMartel and Shapiro! (J1998); Yepesl (2001). 
These variables are remarkable as their use almost completely (completely 
for ideal gas with 7 = 5/3) eliminates explicit dependence on cosmology 
in the model equations. The equations can therefore be integrated using 
standard solvers used in computational fluid dynamics, with no need for 
additional coefficients and corrections. 

The main features of the gasdynamics implementation fo llow closely the 

algorithm used in YK 3 . A second-order Godunov-type solver ( Golella and WoodwardI 
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1984) is used to compute numerical fluxes of gas variables through each cell 



interface, wi th 'left' a n d 'rig ht' states estimated using piecewise linear re- 
construction van Leer! (jl97?l ). For cells that have neighbours at the same 
tree level, the entire integration procedure is formally second order accu- 
rate both in space and in time. For cells that have neighbouring leaves at 
different levels, the accuracy reduces to first order. 

Gasdynamics is coupled to the dynamics of dark matter through the 
common potential, which is used to compute accelerations for both the DM 
particles and the gas. The code employs the same data structures and sim- 
ilar refinement strategy as the FTT used in the TV-body part. However, 
in addition to the DM density criteria it also allows for additional refine- 
ment based on the local density of gas, as well as shock or sharp gradient 
indicators. The refinement criteria can be combined with different weights 
allowing for a flexible refinement strategy that can be tuned to solve the 
Euler equations for an inviscid flow according to the needs of a particular 
simulation. 



2.1.3 Gadget 

In recent years there has been considerable effort in developing tree codes 
that could run in parallel computers. This can alleviate the two most impor- 
tant drawbacks of these algorithms: Large memory requirements and long 
computations per timestep. In the present work we have extensively em- 
ployed the Tree-SPH code Gadget (a rather strange acrony m for Galaxies 
with -D ark matter and Gas int-BracT), recently developed bv lSpringel et al 
and publicly available from this Internet web page 



http : / / www . mpa-garching . mpg . de/gadget / 

Gadget computes gravitational forces wi th the hierarchical octal-tree 
algorithm proposed by Barnes and Hut] l)l98fil ). The collisional fluid (gas) 



is represented by means of Smoothed Particle Hydrodynamics (SPH). The 
code is fully adaptive both in force computation and in time stepping (i.e. 
each particle has its individual mass and time step). Five different particle 
species can be defined by the user, each of them with its own gravitational 
smoothing, which is implemented as a spline-softened potential for point 
masses. Also, there are different criteria for selecting the optimal timestep. 
In what follows we will briefly describe how this code deals with gravity and 
gas dynamics. 



Gravity 

As we said above, gravity forces are computed by means of the BH tree al- 
gorithm. In this scheme, the particles are arranged in a hierarchy of groups 
according to their position in space. When the force on a particular par- 
ticle is computed, the force exerted by distant groups is approximated by 
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their lowest multipole moments. In this way, the computational cost for a 
complete force evaluation can be reduced to order O(NlogN). The forces 
become more accurate if the multipole expansion is carried out to higher or- 
der, but eventually the increasing cost of evaluating higher moments makes 
it more efficient to terminate the multipole expansion and rather use a larger 
number of smaller tree nodes to achieve a desired force accuracy. As a com- 
promise between accuracy and speed, the multipole expansion is terminated 
after the quadrupole moments have been included. 

Force computation then proceeds by walking the tree, and summing up 
appropriate force contributions from tree nodes. In the standard BH tree 
walk, the multipole expansion of a node of size I is used only if 

r > \, (2.9) 

where r is the distance of the point of reference to the centre-of-mass of the 
cell and 6 is a prescribed accuracy parameter. If a node fulfills the criterion 
(|2.9|) . the tree walk along this branch can be terminated, otherwise it is 
'opened', and the walk is continued with all its daughter nodes. For smaller 
values of the opening angle, the forces will in general become more accurate, 
but also more costly to compute. 

The standard BH opening criterion tries to limit the relative error of 
every particle-node interaction by comparing a rough estimate of the size 
of the quadrupole term, ~ Ml 2 /r 4 , with the size of the monopole term, ~ 
M/r 2 . The result is the purely geometrical criterion of equation (|2.9|) . In the 
Gadget code, the acceleration of the previous timestep is used as a handy 
approximate value for the accuracy of force computation. Therefore, it is 
required that the estimated error of an acceptable multipole approximation 
is some small fraction of this a total force. 

The tree-construction can be considered very fast in both cases, because 
the time spent for it is negligible compared to a complete force walk for all 
particles. However, in the individual time integration scheme only a small 
fraction M of all particles may require a force walk at each given timestep, 
and hence the full tree is only reconstructed after a fixed number of timesteps 
has elapsed. 

Time integration of Newton's equations is a variant of the classical Leap- 
frog scheme in which an explicit prediction step is introduced to accommo- 
date individual timesteps for particles. There are different possibilities for 
the choice of timestep size in Gadget. One is based on the second order 
displacement of the kinetic energy, assuming a typical velocity dispersion 
a 2 for all particles. This results is At = a<r/|a|. The other method is 
to constraint the second-order particle displacement (At oc l/y|a|), or a 
more sophisticated one (and more computational expensive) based on an 
estimation of the local dynamical time ( At oc l/\/Gp ) . Ac cording to the 
convergence studies recently carried out bv lPower et al.l l|2002h , we have cho- 
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sen the timestep criterion based on second order particle displacement, with 
the optimal parameters as given in that study. 

Gas dynamics 

Gadget uses an implementation of SPH to estimate pressure forces acting 
on the collisional fluid. SPH is a powerful Lagrangian technique to solve 
hydro dynamical problems with an ease that is unmatched by grid based 
fluid solvers (see iMonaghanl f992 ,for an excellent review). In particular, 
SPH is very well suited for three-dimensional astrophysical problems that 
do not crucially rely on accurately resolved shock fronts because, as we 
mentioned earlier, numerical artificial viscosity introduces some unphysical 
dissipation in the fluid. 

Unlike other numerical approaches for hydrodynamics, the SPH equa- 
tions do not take a unique form. Instead, many formally different versions 
of them can be derived. Furthermore, a large variety of recipes for spe- 
cific implementations of force symmetrisation, determinations of smoothing 
lengths, and artificial viscosity, have been proposed. 

Most of the SPH implementations integrate the standard mass, momen- 
tum and energy conservation equations for the fluid. The computation of 
the hydrodynamic force and the rate of change of internal energy proceeds in 
two phases. In the first phase, new smoothing lengths hi are determined for 
the active particles (these are the ones that need a force update at the cur- 
rent timestep, see below), and for each of them, the neighbouring particles 
inside their respective smoothing radii are found. The Lagrangian nature of 
SPH arises when this number of neighbours is kept approximately constant. 
This is achieved by varying the smoothing length hi of each particle accord- 
ingly, leading to a constant mass resolution independent of the density of 
the flow. 

From the neighbours of each active particle, the density and the rest of 
quantities can be found using an spline interpolation kernel 

N 

f H = J2m j W(T ij ;hi), (2.10) 
j'=i 

where rj,- = r, — Tj. For the rest of the particles, values of for density, 
internal energy, and smoothing length are predicted at the current time 
based on the values of the last update of those particles. An equation of 
state for an ideal gas gives the pressure as a function of density and internal 
energy: Pj = (7 - 1)/^. 

SPH is a slowly converging method with the number of particles. The- 
oretically, SPH equations will be exact in the limit N$ph — ► 00. In cos- 
mological problems, in whi ch large volumes are often simulated, the mass 



resolution is rather coarse. iHernauisti (|1993h has shown that while integra- 



tion of energy equation with SPH formalism results in rather good energy 
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conservation, the entropy is not conserved even for purely adiabatic flows. 
He also showed that the cause of the errors can be found in the use of vari- 
able smoothing and the neglect of relevant terms in the dynamical equations 
associated with the space dependence of hi. 

However, even if fixed smoothing is used in a simulation, simultaneous 
conservation of entropy and energy can only be achieved in the limit of large 
Nsph- Most people believe that violation of entropy is a fair price to pay 
for energy conservation. Therefore, most of SPH co des developed so far do 



show this problem. As have been shown recently ([Springel and Hernquisn 



2002ah . integrating the thermal energy equation with SPH with not enough 



number of particles can lead to significant violations of entropy conservation 
in situations which are important for cosmological simulations of galaxy 
formation. As we will discuss in detail in Section 14.11 this problem with 
entropy dissipation of SPH codes could be the reason o f the discrepancie s 
that showed up in the Santa Barbara cluster comparison (|Frenk et aljf l999^. 



Two possibilities have been proposed to overcome this problem: One is to 
include in the formulation of the dynamical equations for SPH particles the 
V/i terms associated with spatial derivatives of the smoothing length (e.g. 
the AP 3 M-SPH code DEVA l Serna et all2002h . The other solution, proposed 
by Springel and HernquistJ lj2002al ). is an elegant formulation of SPH which 



employs variable smoothing lengths and explicitly conserve both energy and 
entropy. To this end, the equations of motion for the fluid particles were 
derived from the Lagrangian defined as: 

N N 

£(q, q) = 2 E m ^ ~ — E mtAipT 1 (2-H) 
i=i ~ i=i 

in the independent variables q = (ri, . . . , rjv, h\, . . . , hjy), where the thermal 
energy acts as the potential generating the motion of SPH particles. The 
densities pi are functions of q defined by the standard SPH kernel inter- 
polation from particle positions. The internal energy is computed in terms 
of the specific entropy of fluid elements, A(s), which for an ideal gas is 
u i = ^Ti/°7 _1 - The quantities A{ are treated as constants (i.e. the flow is 
assumed to be strictly adiabatic). 

The smoothing lengths hi are selected by requiring that a fixed mass is 
contained within a smoothing volume, which provides N constraints 

4"7T o 

&(q) = —h\ Pi - M sph = (2.12) 

on the coordinates of the Lagrangian (note that, for SPH particles of the 
same mass, the requirement of constant M sp h is equivalent to a fixed number 
of neighbours). The equations of motion reduce thus to 



A N 

dvi 



dt 



fJ&iWaihi) + fj^ViWijihj 



(2.13) 
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I I K dpi 
' 3pi dhi 



and the abbreviated form 



where the fa are defined by /, 

Wij(h) = W(\ri — Tj\, h) has been used. 

Because the potential (thermal) energy of the Lagrangian depends only 
on coordinate differences, the pairwise force in equation Q2.13[l is automat- 
ically anti-symmetric. Total energy, entropy, momentum, and angular mo- 
mentum are therefore all manifestly conserved, provided that the smoothing 
lengths are adjusted locally to ensure constant mass resolution as defined 
by equation (|2.12j) . Moreover, the V/i terms are implicitly included in the 
equations through the factors. 

As in the standard (energy) implementation of SPH, artificial viscosity 
must be included to allow for the handling of shocks. A viscous force 



A N 

dvj 
dt 



^mjUijViWy , (2.14) 



is added to the acceleration given by equation ()2.13j) . The resulting dissi- 
pation of kinetic energy is exactly balanced by a corresponding increase in 
thermal energy if the entropy is evolved according to 

^ = --^y^CipijUi) + ^^pr^mjllijVij ■ ViWij, (2.15) 



Pi 



=i 



which assumes that entropy is only generated by the artificial viscosity in 
shocks, and by possible external sources of heat ((C(pi,Ui)), if they are 
present (e.g. feedback from stars). 

The modifications needed in the original Gadget code to account for 
this new version of SPH are not very important. A slight alteration is 
required for the algorithm which updates smoothing lengths, since it is now 
necessary to ensure that a 'constant mass' resides within a smoothing volume 
rather than a constant number of neighbours. 

Tests performed by ISpringel and Hernauistl fe002al ) with this new im- 



plementation of SPH show that entropy is much better conserved, avoiding 
the substantial overcooling produced in poorly resolved haloes. In the stan- 
dard SPH implementation based on the thermal energy equation, the lack 
of resolution can lead to insufficient compressional heating in the accretion 
flow, resulting in severe entropy losses. 

For our numerical experiments of cluster formation we have used the 
new version of Gadget , based on the entropy-conserving implementation 
of SPH. Nevertheless, in order to check the differences with respect to the 
standard SPH technique, we have run several test simulations of cluster 
formation with both versions of the parallel code: The standard, public 
available, version and the entropy implementation of SPH, kindly provided 
by Volker Springel. Moreover, to compare SPH results with other numerical 
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integration schemes, some simulations of the same initial conditions have 
been rerun with the Eulerian gasdynamical code ART by Andrey Kravtsov 
( Kravtsov et al. 20021 : Nagai and Kravtsovll2002l ) . The outcome of this com- 
parison is reported in Section f4. II 



2.2 Analysis 

The output of a numerical experiment is a data file with the positions and 
velocities of dark matter and gas particles, gas densities and temperatures, 
and so on. In order to make a physical interpretation of these data, the 
analysis procedure is as important as the numerical simulation itself. 

In this section, we describe the conventions followed in the analysis of 
the results of our numerical simulations. We address the questions of how 
to set the resolution limit, the definition of objects from particle data, and 
the computation of several magnitudes of physical interest, as well as their 
radial profiles. Throughout the present work, all magnitudes are expressed 
in units of Hq = lOO/i km s _1 Mpc -1 . 



2.2.1 Resolution limit 



Maybe the most common problem that must be faced by numerical simula- 
tions is the lack of mass and force resolution. Spurious effects arising from 
a low number of particles constitute a major source of uncertainty in any 
numerical experiment. In galaxy cl usters, poor resolution leads to the so - 
called overmerging problem (see e.g. Iwhite et al.lll987 : Klvpin et al. 199fll ). 
namely the loss of substructure within the dark matter halo. Moreover, res- 
olution is of critical importance in determining the radial profiles of both 
dark matter and baryons at very small radii. 

It is not clear which numerical effects may determine the minimum scale 
above which the results of a given code can be trusted, but it is quite likely 
that this convergence scale is determined by a complex interplay of all pos- 
sible numerical effects. Although there have been some recent attempts at 
unravelling the role of numerical parameters on the structure of simulated 
dark matter haloes (e.g. Moore et al . 199 8al:lKnebe et alJ koOO: Klvpin et al 



hflflflL I2OO1I : lOhiema, et al.1 1200(1 IPower et al J 1200^ . the conclusions fr 



om 



these works are still pr eliminary and, in so me cases, even contradictory. 

To cite an example. iMoore et al. ( 1998ah argue that the smallest resolved 
scales correspond to about half the mean inter-particle separation within the 
virial radius, and conclude that many thousands of p articles are ne e ded t o 



resolve the inner density profile of dark matter haloes. iKlvpin et al. (2001), 



on the other hand, conclude that the convergence scale is ~ 4 times the 
formal force resolutio ns, or the rad i us con taining 200 dark matter particles, 



whichever is larger. iGhigna et al.l (J2000) suggest a similar criterion based 
on the gravitational softening length, and argue that convergence is only 
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achieved on scales that contain many particles and that are larger than 
about ~ 3 times the scale where pairwise forces become Newtonian. 

Understanding the precise role of numerical parameters is clearly needed 
before a firm theoretical prediction for the structure of CDM haloes on ~kpc 
scales may emerge. This question is particularly difficult because of the lack 
of a theory with which the true structure of dark haloes may be predicted 
analytically, so the best that can be done is to establish the conditions under 
which the structure of a simulated dark halo is independent of numerical 
parameters. 

Furthermore, we note that most convergence criteria have been derived 
considering the minimum scale at which the density profiles of different 
resolution simulations are consistent among themselves. Other properties 
of dark matter haloes, such as their three-dimensional shape, their detailed 
orbital structure, or the mass function of substructure haloes, may require 
different resolution limits. For the radial profiles related to gasdynamical 
quan tities, the convergen ce scale has been investigated in several studies 



(e.g. Borgani et al. 20021 '). The numbe r of SPH neighbou rs plays a crucial 
role in the accuracy of results (see e.g. IPower et al.l 12002). In all our SPH 
simulations, this parameter has been set to 40 particles. 

Having all these considerations into account, we adopted the follow- 
ing criteria to set the resolution limit of the numerical experiments using 
particle-based codes (Gadget) described in this chapter: 

• 100 gas particles contained within r m ; n 

• 200 dark matter particles 

• r min > 4e, where e is the gravitational softening length 1 

Usually, the first condition imposes the most stringent constraints on 
our resolution, since gas tends to be substantially less concentrated than the 
dark matter due to the repulsive pressure force. Hence, the radial profiles 
computed for Gadget feature an apparent lower spatial resolution than 
those of the N-body code ART, but this only reflects the presence of gas in 
the former simulations. 

This resolution limit is only applied to the minimum scale that has been 
considered reliable in the radial profiles. Global quantities (see Section 12.2.41 
below) have been computed from all the particles within the virial radius of 
the clusters. 

2.2.2 Object finding 

Obviously, the first step that must be taken to analyse a cluster of galaxies 
(or any other object) is to locate it. Identification of haloes in high density 

1 In Gadget, using a spline-based smoothing potential for point particles, pairwise 
force becomes Newtonian for r ~ 2.3e 
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environments, such as groups and clusters of galaxies, is a challenging prob- 
lem. Currently, there are several methods to accomplish this task, each one 
with its own merits and drawbacks. To mention some o f the most popular, 
the Friends Of Friends (FOF, see e.g. iDavis et alJll98.4 ) algorithm connects 
those particles that are separated less than a given linking length. A h ier- 



archical version of this method has been proposed ( Klvpin et al. 1999h to 



overcome the merging of apparently distinct haloes when the linking length 
is set too large and the missing of many particles w hen it is too small. An 
alternative is the ' spherical overdensity' scheme (e.g. Lacev and Cole! 19941 : 
Klvpin et al. 1996^ . which looks for local maxim a of the density field a bove 



some threshold defined as the 'virial' overdensity. ISpringel et al.l ( 2001al ) use 
a more elaborate algorithm that limits the objects by the isodensity contour 
that traverses a saddle point. 

For our analysis, we have chosen the Bound Density Maxima (BDM, see 
e.g. lCohn etaD ll999: Klv pin et al.lll999l ^ galaxy finding algorithm to iden- 
tify bound structures from particle positions and velocities. This approach 
is based on the 'spherical density' outlined above, but particles that are not 
bound to the identified halo are removed recursively until a self-consistent 
definition of the object is found. The BDM algorithm proves extremely help- 
ful in finding substructure within the dark matter haloes of galaxy clusters. 
As we discuss in Section \'2. 3. ll below. the presence of significant substructure 
within the virialised region of a cluster is used as an indicator of dynamical 
activity. 

Another important issue is the determination of the centre of mass. This 
question is particularly relevant in those cases where the cluster halo is not 
exactly spherically symmetric, or when the study focuses on the detailed 
properties of the innermost regions of the simulated clusters. 

We use as a first approximation the position given by the BDM algo- 
rithm, and then an iterative technique is employed in which the centre of 
mass is found for particles contained within a shrinking sphere. At each 
step of the iteration, the centre of the sphere is reset to the last computed 
barycenter (as a starting point, we use the BDM value) and its radius is 
reduced by 5%. The centre of mass is thus computed recursively until a 
convergence criterion is met. In our case, the iteration is stopped at the res- 
olution limit described above, checking that the variation in the final centre 
of mass is less than 1 h" 1 kpc. 

Halo centres identified with this procedure are quite independent of the 
parameters chosen to start the iteration, provided that the initial sphere 
is large enough to encompass a large fraction of the system. In a multi- 
component system, such as a dark halo with substructure, this algorithm 
isolates the densest region within the largest subcomponent. In more regular 
systems, the centre so obtained is in good agreement with centres obtained 
by weighing the centre of mass by the local density or gravitational potential 
of each particle. 
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2.2.3 Radial profiles 

Throughout most part of this work, our analysis of the physical properties of 
galaxy clusters will focus on the spherically-averaged radial profiles of several 
quantities of physical interest. Although the detailed structure of clusters 
displays significant departures from spherical symmetry, radial profiles offer 
a simple description of the dynamical state of the system, providing a great 
deal of helpful information about both the dark and baryonic components. 

This section is devoted to clarify those aspects of our analysis that are not 
straightforward. The spherically-averaged profiles are measured by binning 
the particles according to their distance from the centre of mass of the at 
z = 0. Logarithmic bins in steps of 0.05 dex are used in all cases, except 
in the computation of the radial profiles of projected quantities, for which 
a constant step of 10 h^ 1 kpc is used. In this case, the cluster properties 
are computed within cylindrical shells oriented along the coordinate axes, 
truncated at a distance of 3 h^ 1 Mpc. The final profile is taken to be the 
mean value of the three projections in order to reduce noise. 

Dark matter 

Density, cumulative mass and circular velocity (V^ 2 (r) = GM(r)/r) are com- 
puted the usual way. In order to study the kinetical structure of our clusters 
of galaxies, we have obtained several radial profiles related to the velocity of 
the particle with respect to the centre of mass. For each bin, the movement 
have been decomposed into the radial and tangential directions, as well as 
in bulk (i.e. average) and random motions (of (r) = (vf) — {vi) 2 ). Angular 
momentum has been computed from the tangential component of the ve- 
locity, assuming that the radius is constant for all the particles in the bin. 
Ordered rotation is computed from (vt g ), while the contribution of random 
motions comes from the tangential velocity dispersion. 

The orbits of individual particles have also been classified according to 
their apocentric radii (see Section 13.2.2)) . This is done by computing the 
gravitational potential (j)(r) — (po = J Q GM(x)/x 2 dx, assuming spherical 
symmetry. This approximation incurs in some uncertainty due to the time 
variation in the mass distribution of the cluster. This effect is more im- 
portant in the outer parts, where the cluster is still accreting a significant 
amount of matter at the present epoch. 

Once the gravitational potential is known, the pericentric and apocentric 
radii that a particle reaches during its orbit are easily computed from its 
position and velocity. A useful quantity is the eccentricity, defined as e = 
(r max — r m i n )/(r max + r m i n ). It measures the depth that particles are able 
to penetrate into the cluster potential well; a value e = corresponds to 
perfectly circular orbits, whereas e = 1 indicates that the particle trajectory 
goes through the centre. 
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Gas 

Gas density and mass are computed in an analogous manner as for the 
dark matter. We will use the notation Ff, to denote the global baryon 
fraction, Fb = M g /M. The local baryon fraction will be referred to as 
fb(r) = Pg{r)/ pdm( r )- Note that the total mass is used in the definition of 
Ft, while fb is the ratio of the baryonic to the dark matter density. The 
difference, though, is quite small. 

Temperature of gas particles is inferred from its internal energy. Bolo- 
metric X-ray luminosity is approximated by pure bremsstrahlung emission, 
where the emissivity is given by rigA(T), that is, the square of the electron 
density n e = p g /(pm p ) times the cooling function A(T) (m p is the proton 
mass and a = 0.65 is the mean molecular weight of a fully ionised plasma). 
Following Na varro et al we will assume that, at the temperatures 



relevant for galaxy clusters, the cooling function can be well approximated 
by 

/ T \ 1/2 

A(T) ~ 1.2 x 10~ 24 erg s _1 cm 3 (2.16) 



,1 keVy 

and hence the luminosity of a set of SPH particles can be computed as 

N 

-24 \^ m iPi rpl/2 



L X = 1.2 x 10- 24 V -p%Tl' 2 (2.17) 
f^(pm p y 



where mi, pi and Tj are the mass, density and temperature of each gas 
particle. 

Other quantities can be obtained from the radial density and temper- 

diog(p) 

entropy, defined as S{r) = Tn e 



-2/3 



Average profiles 

A convenient way to represent our results is to combine the information 
obtained from a certain number of clusters into a single plot. Our approach 
is based on the following scheme: 

1. Compute individual profiles 

2. Define radial bins (points to plot) 

3. For each bin, interpolate the value of every profile at that point 

4. Compute mean value and standard deviation 

Averaged profiles reported in the present work have been computed ac- 
cording to this simple prescription. The quoted error bars correspond to the 
standard deviation (note that, in some cases, this quantity can be dominated 
by Poisson noise due to the low number of clusters considered). 
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2.2.4 Global properties 

In many occasions, we will be interested not only in the radial structure of 
galaxy clusters, but also in their bulk properties. Therefore, we must define 
their average value up to some radius that we define as the cluster limit. 

As it happens in observations, there is no unambiguous prescription to 
define this limit. In a Standard Cold Dark Matter (SCDM) cosmology, 
spherical collapse theory predicts that dark matter haloes should be in virial 
equilibrium at densities of roughly 200p c . The radius at which the cumu- 
lative overdensity reaches 200 has often been used to mark the boundary 
of both observed and simulated objects. However, this estimate is a very 
simple approximation to the non-linear process of virial relaxation. As we 
will show in Section I4.2.H the assumption that clusters of galaxies are in 
equilibrium at these large radii holds only marginally. 

In a ACDM universe, this simple prescription based on spherical col- 
lapse yields an even lower density for virialisation (~ 100p c ). However, for 
historical reasons, it is not uncommon that the value corresponding to the 
SCDM model is used to define the limit of clusters. Since the difference in 
most physical quantities computed at noo or r2oo is n °t large, and the latter 
is most often quoted in observational studies, throughout this thesis we will 
use the term 'virial' to denote quantities defined at A = 200, that is, 

M 200 = y200p c r| 00 ( 2 -!8) 

In any case, the overdensity at which a magnitude is evaluated will always 
be indicated by a numerical subscript. For instance, in Section 14.31 the 
scaling relations between the total mass, luminosity and temperature are 
investigated at several values of the overdensity threshold A, similar to those 
observed by the current generation of X-ray satellites. 

In order to make a more consistent comparison with these observations, 
we compute the emission-weighted temperature of our galaxy clusters as 

sr^N T 3/2 
2_/i=l Pi ± i 

Unless stated otherwise, 'temperature' or T' will always refer to the 
mass- weighted temperature, T = ^2iL\Ti/N. All other quantities are com- 
puted in the obvious way. 

2.3 Description of the simulations 

Most of the numerical analyses reported in this work is based on the study 
of a sample of 15 galaxy clusters simulated both with the adaptive Eulerian 
code ART (including only dark matter particles) and the e xplicit entropy- 



conserving implementation of the Lagrangian code Gadget (jSpringel and Hernquist 
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2002ah . which uses a variant of the SPH technique to account for the gas- 



dynamics of the baryonic component of the intracluster medium. 

In addition to this cluster sample, we have considered also th e initial con- 



dition s used in the Santa Barbara Cluster Comparison Project (|Frenk et al 



1999), as well as a different set of simulations (accomplished with the Eu- 



lerian YK 3 ) that allow us to investigate the star formation and feedback 
effects in cluster and group environments. 



2.3.1 Cluster sample 

Our sample comprises a total number of 15 clusters of galaxies, selected 
from a low-resolution simulation with 128 3 dark matter particles. In total, 
we have performed 7 independent numerical experiments running Gadget 
on a SGI Origin 3800 parallel computer at Ciemat (Spain), using 32 CPU 
simultaneously. The average computing time needed to run each simulation 
was ~ 8 CPU days (6 x 10 5 s). The same clusters have also been simulated 
with the N-body version of ART on the Hitachi SVR at the LRZ (Germany). 

The properties of the cluster sample are thoroughly discussed in Sec- 
tion !3.21 devoted to the dark matter distribution within the simulated haloes. 
In Chapter 01 the structure and scaling relations of the ICM gas for these 
clusters are investigated. 



Initial conditions 

In a cubic volume of 80 Mpc on a side, an unconstrained realisation 
of the power spectrum of density fluctuations corresponding to the most 
favoured ACDM model (Q m = 0.3, A = 0.7, h = 0.7 and a 8 = 0.9) 
was generated for a total of 1024 3 Fourier modes. The density field was 
then re-sampled to a grid of 128 3 particles which were displaced from their 
Lagrangian positions according to Zeldovich approximation up to z = 49. 
From this initial conditions, ART (dark matter only) and Gadget (dark 
and gas) codes were used to evolve the particles up to z = 0. 

Selected clusters have been re-si mulated with higher resolution by means 
of the multiple mass technique (see Klvpin et alJ 200 ll . for further details). 



For each one of them, we have computed the particles in a spherical re- 
gion around the centre of mass of the 128 3 low-resolution counterpart at 
z = 0. Mass resolution is then increased by adding smaller particles in 
the Lagrangian volume depicted by these particles, including the additional 
small-scale waves from the ACDM power spectrum in the new initial condi- 
tions. 

This is illustrated in Figure l2~Tl where the colour code is used to indicate 
the refinement level. A buffer of medium-resolution particles (blue) avoids 
undesired two-body encounters between high-resolution (red) particles and 
the coarsest level (black) of refinement. 
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Figure 2.1: Mass refinement technique 



Some of the clusters in the sample described in the present section were so 
close to each other that their Lagrangian spaces overlapped at high z. Hence, 
these objects have been simulated together and only 7 experiments have 
been necessary in order to re-simulate the full sample with high resolution. 
Individual images of every cluster can be found in Appendix 1X1 

Parameters of these simulations are summarised in Table l2~T1 where the 
first column indicates the clusters that where refined in that run. In the 
second column, the total number of particles in the 80 Mpc box is 
given, and the third column states the number of high resolution particles 
for each type of matter (i.e. N/^ dark matter particles and N/w gas particles 



Clusters 


N 


N hi 


JU DM 


M hi 

-"-•gas 


A 


4504797 


1190016 


29.6 


2.12 


B 


4719315 


1291008 


29.6 


2.12 


C 


5454182 


1651648 


29.8 


2.01 


D, H 


4411785 


1136064 


29.6 


2.12 


E, L, M 


4941262 


1389440 


30.0 


1.80 


F, G, K 1; K 2 


4585734 


1217088 


29.6 


2.12 


Ij Ji, J2 


5032957 


1438848 


29.6 


2.12 



Table 2.1: Total number of particles in each simulation, number of high resolution 
particles (Nhi dark matter + N/jj gas) and their mass in 10 7 /i _1 M . 
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Cluster State r 200 M 2u0 L% 00 T 200 T^ 00 

A Minor 931 18.96 70.06 2.000 2.873 

B Minor 871 15.57 20.30 1.860 2.620 

C Relaxed 871 15.53 42.65 1.958 2.810 

D Minor 771 10.77 11.49 1.301 1.642 

E Relaxed 719 8.74 12.53 1.287 1.866 

F Major 661 6.79 12.47 1.059 1.195 

G Major 638 6.10 4.50 0.844 0.818 

H Relaxed 618 5.56 16.60 1.107 1.614 

I Major 581 4.60 2.34 0.666 0.666 

Ji Relaxed 584 4.67 9.99 0.937 1.498 

Ki Relaxed 557 4.06 4.98 0.690 1.043 

L Minor 547 3.84 1.21 0.624 0.779 

M Major 503 2.99 0.64 0.545 0.610 

K 2 Minor 497 2.89 3.08 0.470 0.808 

J 2 Relaxed 491 2.77 8.22 0.673 0.979 



Table 2.2: Physical properties of the clusters at z = 0. r2oo in h~~ Y kpc, enclosed 
mass in lO 13 /)." 1 M Q , bolometric X-ray luminosity in 10 25 h erg s _1 and 
average temperatures (mass and emission-weighted) in keV. One or two 
asterisks beside the cluster name indicate minor or major merging activity 
(see text). 



have been used in each experiment). The last two columns show the mass 
resolution in units of 10 7 ft, _1 Mq, corresponding to the three levels of mass 
refinement that have been employed in all cases. 

The gravitational softening length has been set to e = 2 — 5 h~ l kpc, 
depending on the number of par ticles within the virial radius of each indi- 
vidual cluster ( Power et al.ll2002l ) . The minimum smoothing length for SPH 
was fixed to the same value as e. 



Physical properties 

Table displays the physical properties of the simulated clusters at z = 0. 
As discussed in Section 12.2. 4[ the subscript '200' refers to the over density 
with respect to the critical value p c ~ 2.8 x 10 n /i 2 M Mpc -3 (not to be 
confused with the mean density of the universe, p m = Q m p c ). 

Clusters have been sorted (and named) in Table l2~2*l according to their 
virial mass at the present day. Clusters J 2 and K 2 are an exception to this 
rule, since they are two small groups falling into Ji and Ki respectively. As 
a matter of fact, the smallest objects in the present sample could as well 
be considered rich groups instead of poor clusters. Not existing a clear-cut 
distinction between these categories, we will apply the term 'clusters' to all 
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the objects regardless of their virial mass. 

In order to study the effect of mergers and close encounters, we labelled 
as minor merger any cluster in which we are able to identify a companion 
structure within r2oo whose mass is greater than 0.1 M200; if the mass of 
the companion is above 0.5 M200, then we classify the merger as major, 
otherwise, the cluster is assumed to be a relaxed system in virial equilibrium. 

The results of this classification scheme are quoted in the second column 
of Table E21 Clusters named Jj and Kj are relatively close pairs, but they 
are separate enough (~ 2 — 3 Mpc) not to be considered as mergers. In any 
case, it is remarkable the amount of dynamical activity found in our sample, 
specially for low-mass systems. This is in agreement with the results of 
Oottlober et al. (|200ll ). who found that the typical merging rate in groups 



is much higher than in clusters, even for galaxies of the same mass. 

Another interesting issue is that minor mergers tend not to show signifi- 
cant features in the X-ray images, although some of their physical properties 
can be noticeably affected (see Appendix ^J. High-resolution observation 
(such as those based on XMM or Chandra satellites) are required in order to 
detect asymmetries in the inner parts of most of our merging clusters. The 
consequences concerning observational results based on the assumption of 
hydrostatic equilibrium are difficult to evaluate. 



2.3.2 Santa Barbara Cluster 



The Santa Barbara Cluster Comparison Project ( Frenk et al. 19991 ) was a 



compilation of the results of different hydro dynamical codes, applied to the 
formation of an X-ray cluster of galaxies in a SCDM universe. In order 
to test the two different formulations of SPH implemented in Gadget, we 
have used both versions of the code to carry out simulations with these 
initial conditions. 

Numerical experiments described in this section have been carried out 
in the SGI Origin 2000 at the Centro Europeo de Paralelismo de Barcelona 
(CEPBA), the Astrophysicalisches Institut Potsdam (AIP) and the Na- 
tional Center for Supercomputing Applications (NCSA). Our findings are 
presented in Section E~T1 



Initial conditions 

The initial conditions for the Santa Barbara cluster experiments have been 
supplied by the authors, both in Eulerian and Lagrangian format. They are 
publicly available on the Internet, at 

http : / / star-www. dur . ac.uk/~csf/clusdata/ 

The initial fluctuation spectrum was taken to have an asymptotic spec- 
tral index n = 1, and shape parameter T = 0.25. The cosmological pa- 
rameters are those of the SCDM model, that is Q m = 1, = 0, h = 0.5 
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Version 


iVeff 


MM 


e [/i" 1 


Hj 


04 


1NO 


on 


T? 

rL 


1 Oq3 


1NO 




T71 

E 


128 


AT 

JNo 


5 


E 


128 3 


Yes 


15 


S 


128 3 


No 


20 


S 


128 3 


No 


5 


S 


128 3 


Yes 


15 


S 


256 3 


Yes 


2 



Table 2.3: Santa Barbara Cluster simulations 

and f^b = 0.1. The cluster perturbation was chosen to correspond to a 3a 
peak of the primordial density field smoothed with a Gaussian filter of ra- 
dius rn = 10 Mpc (in exp[— 0.5(r/rrj) 2 ]). The perturbation was centred on a 
cubic region of side L = 64 Mpc. 

This cluster has been simulated several times, varying the formulation of 
SPH used in the integration of the gasdynamic equations, the total number 
of particles, and the gravitational softening length. 

Table |2~31 summarises the values of these parameters. The implementa- 
tion of SPH is quoted by the letter in the first column, where 'E' stands for 
'energy', and 'S' for 'entropy'. In each case, the corresponding conservation 
equation is used to evolve the gas temperature in time. The second column 
quotes the effective resolution of the simulation, and the third, whether 
multi-mass techniques have been applied or not. When multiple masses 
have been used, the effective resolution correspond to the highest level of 
refinement; otherwise, N e g equals the total number of particles in the sim- 
ulated box (actually, one half, because there are as many gaseous as dark 
matter particles). The gravitational softening length is quoted in the last 
column. 



Physical properties 

Quite expectedly, given that all the simulations correspond to the same clus- 
ter, the final global properties of the system are similar in all the experiments 
quoted in Table fOI r 200 - 1-35 h~ l Mpc, M 200 ^ 0.58 x 10 15 h~ l Mpc, 
T 20 o = 4.6 keV, Lf 00 ~ 1.27 x 10 45 h erg s" 1 



Th ese figures are consistent with the values reported in iFrenk et al 
(1999) for the average of all the numerical experiments presented in the 
Santa Barbara study. However, the detailed structure of the radial profiles 
of the diffuse gas component shows systematic trends as a function of the 
SPH implementation. A detailed discussion of this effect, as well as that of 
poor resolution, will be given in Section EU 
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ACDM 0.325 0.025 0.65 0.7 0.9 
SCDM 0.95 0.05 0.5 1.2 
BSI 0.95 0.05 0.5 0.6 

Table 2.4: Cosmological models 

2.3.3 Star-forming clusters 

In order to investigate the formation of stars and the effects of the environ- 
ment on the star formation rate (SFR) of the simulated galaxies, we have 
performed a series of numerical experiments with the code YK 3 , described 
in Section [2.1.11 The numerical code has been parallelised using OpenMP 
compiler directives. Therefore, it runs very efficiently in parallel computers 
with either Shared or Non Uniform Memory Access (NUMA) architectures. 
The simulations have been performed on a variety of machines. 

These simulations were aimed in principle to a broader goal than the 
study of galaxy clusters, namely the evolution of the SFR density in different 
cosmologies and environments. Therefore, they include many small-volume 
numerical experiments in which no clusters of galaxies are formed. These 
experiments, however, offer interesting insights on the star formation history 
of the universe, as well as the dramatical changes that the star formation 
activity of a galaxy undergoes when it is accreted into a massive group or 
cluster of galaxies. These issues are covered in detail in Chapter 

Initial conditions 

As has been discussed in Section II. If the ACDM model has proven to be 
very successful in describing most of the observational data at both low 
and high redshifts. However, in order to study the influence of cosmology 
on the global SFR, we have investigated two other scenarios: the standard 
CDM model (dominated by dark matter density), and the BSI model with 
a Broken Sc ale Invariant perturb ation spectrum (as predicted by double 
inflation, see iGottldber et alJll99lh . The main parameters describing these 



cosmologies are summarised in Table l2~H 

The BSI model assumes a dark matter-dominated universe similar to 
the SCDM cosmology, but primordial fluctuations on small scales are signif- 
icantly suppressed from its power spectrum (the overall normalisation, as in 
the other two models, is set by COBE data on large scales). The resulting 
P(k) is almost the same spectrum as in the rCDM model, in which the pri- 
mordial perturbation spectru m has been changed d ue to the hypothetical 



decay of massive r neutrinos (E fstathiou et al 



ged d m 

Eaa3). 



All the experiments described in this section are based on completely in- 
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Experiment n T N L^ ox L ce a Mdm A V 



ACDM1 


1 


350 3 


30 


85.7 


28.1 


200 


100 


ACDM2 


1 


270 3 


21.4 


79.3 


22.2 


200 


100 


ACDM3 


1 


300 3 


12 


40 


2.85 


200 


100 


ACDM4 


9 


128 3 


5 


39.1 


2.22 


200 


100,500,1000 


ACDM5 


18 


128 3 


5 


39.1 


2.65 


0,50,200 


100 


SCDM 


18 


128 3 


5 


39.1 


3.96 


0,50,200 


100 


BSI 


18 


128 3 


5 


39.1 


3.96 


0,50,200 


100 



Table 2.5: Resolution and feedback parameters. Number of realisations, total number 
of particles/cells, box length in hr 1 Mpc, cell length in h" 1 kpc, mass 
of dark matter particles in 10 6 M Q , supernova feedback parameter and 
overdensity threshold for star formation. 



dependent initial conditions (i.e. they do not belong to a consistent sample, 
as in Section |2.3.1|) . As will be discussed in more detail in Sections 15.1.21 
and 15 .21 the effects of cosmic variance due to a small simulated (or observed) 
volume are not negligible on the estimation of the mean star formation rate 
density. 

Physical properties 

Table l2~5l summaries the set of experiments performed with YK 3 . In order 
to improve statistics, several realisations have been run when the number 
of particles was low enough (N = 128 3 ), changing the random seeds used 
to generate the initial conditions. The number of realisations is given by 
n T in the second column. Next 4 columns show, respectively, the number 
of particles, box size, cell size (which determines spatial resolution in these 
simulations) and mass of each dark matter particle (mass resolution) . In the 
last two columns, we quote the values of the supernova feedback parameter 
A and the overdensity threshold T> required for star formation (details on 
the role of these parameters in YK 3 can be found in Section [2.1.1)1 . 

The comparison of experiments ACDM5, SCDM and BSI allows us to 
determine the effect of cosmology and feedback from supernovae on the 
global SFR history. The photoionisation prescription can be tested basing 
on the results of experiment ACDM4, and the dependence on resolution and 
volume will be thoroughly studied from the whole set of ACDM experiments. 

Finally, experiments ACDM 1 and 2 simulate a volume large enough to 
contain a cluster of galaxies. We will use their results to investigated the 
star formation history of clusters in Section 15.21 focusing on the alterations 
of the star formation activity of galaxies when they fall into the potential 
well of the cluster. 



Chapter 3 

Dark matter 



If the sun refuse to shine 
I don't mind 
I don't mind 

- Jimi Hendrix : If 6 was 9 (1968) - 



^tcoetiing to the current cosmological paradigm, most mass in the uni- 
verse must be in the form of non-baryonic cold dark matter particles, whose 
nature is yet an unknown that has been puzzling both astronomers and par- 
ticle physicists during the past few decades. As far as we are concerned, 
cold dark matter can be considered as a perfect fluid with negligible pres- 
sure, thus obeying the collisionless Boltzmann equation. Not being coupled 
to the electromagnetic field, CDM particles do not emit or absorb light, and 
hence (neglecting some form of self-interaction) gravity is the only physical 
process that must be accounted for. 

In clusters of g alaxies, the baryon frac tion is expected to be close to the 
cosmic value (e.g. White and Frenkl 1991 '). which implies that the amount 
of cold dark matter is higher by approximately one order of magnitude than 
the mass in baryonic form. Structure formation will be thus driven almost 
entirely by the CDM component. Galaxies and clusters are assumed to 
form by cooling and condensation of baryons into the potential wells created 
by dark matter haloes, which were able to collapse during the radiation- 
dominated epoch. 

In this chapter, we will study the gravitational and kinetic structure 
of clusters of galaxies. The radial distribution of mass, density, velocity 
dispersion and angular momentum will be investigated in detail, focusing 
on the existence of common patterns usually known in the literature as 
'universal profiles'. We will try to understand the physical origin of these 
patterns in the dark matter distribution by means of a simple analytical 
treatment based on the secondary infall model. 
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Most of the work refers to the cluster sample described in Section T2.3.1I 
Results concerning the densit y profiles of these haloes and the dependence on 



environment can be found in Ascasibar et al. ( 2002ci ). The orbital structure 



of galaxy clusters, as well as t he numerical treatment o f spherical collapse 



are to going to be published in Ascasibar et al. ( 2002bh . 



3.1 The case for a universal density profile 



One of the most interesting properties of the dark matter haloes found in 
numerical simulations is the apparent universality of their radial density 
profiles, valid over several orders of magnitude in mass. 



Contrary to the early numerical work o flQuinn et all ([19861 ) and lFrenk et al 



(1988), in wh ich dark matter haloes s howed an i sothermal st r ucture (i.e. 
p(r) oc r" " 2 ), iDubinski and Carlberd (H) and ICrone et all JSJ had 
enough resolution to detect the f irst evidence of non-power- law density pro- 
files. Building upon these results, Navarro et all ( 19961 . 1997 , hereafter NFW) 
found that haloes in their numerical simulations could be fitted by a simple 
analytical function with only two parameters (e.g. a characteristic density 
and radius). 

The main result of NFW experiments was that the radial density profile 
was steeper than isothermal (p oc r~ 3 ) for large radii, and shallower (but 
still diverging as p oc r _1 ) near the centre. The corresponding logarithmic 
slope 

dlog p(r) 



a(r) 



dlogr 



(3.1) 



varied smoothly between these two extreme values, being equal to isothermal 
(a = — 2) at the characteristic radius only. 

iNavarro et al. (1997) further show that the characteristic parameters of 
a dark matter halo seem to correlate to a certain extent. Should this be true, 
the final mass distribution of objects of different mass could be described in 
terms of a one-parameter family of analytical profiles. This also implies that 
many relevant quantities, such as virial mass, central density, or formation 
time, must be correlated. 

Similar results have been found in several independent studies, using 
much higher mass and force resolution than the original NFW paper. How- 
ever, there is still some disagreement about the inner most value of t he log- 
arithm i c slope a and its depend ence on resolution. TMoore et all (Il998 
ll999bh : lGhigna et all (|l998l . l200Ch and iFukushige and Makind l|l997l .l: 
find eve n steeper density profiles near the centre (a ~ —1.5), whereas other 



authors ( Jing and Sutolboool : Klvpin et al .11200 ih claim that the actual value 
of a may depend on halo mass, merger history and substructure. 
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3.1.1 Models 



In order to investigate the existence and shape of a 'universal' density pro- 
file, we will compare the mass distribution of our simulated clust ers with the 
functional forms advocated by NFW and IMoore et al. (1999b). Although 
these are not the only profiles proposed in the literature 1 , they are by far the 
most widely referenced. Below, we describe the main features of both pro- 
files, as well as some derived quantities that we will use in the comparison. 
Further details concerning these functions, such as the associated gravita- 
tional potential or t he radial velocity dispersion can be easi ly computed 
analytically (see e.g. Lokas and Mamor]l200ll : ISuto et al.lll998h . 



NFW 

The radial density profile proposed by NFW 

Ps 



p{r) 



,1 + - 



(3.2) 



depends on the characteristic density p s and the scale radius r s of the dark 
matter halo. This profile is often written in terms of the concentration 
parameter 

c 



^200 



as 



p(r) 



200p e c 2 g(c)/3 
s(l + cs) 2 



where we have used the definitions 



s = 



^200 



; 9(c) 



n -l 



ln(l + c) 



1 + c 



(3.3) 
(3.4) 

(3.5) 



and r2oo follows the convention given in Section 12.2.41 (i.e. 200 times the 
critical density of the universe). 

This profile steepens monotonically with radius. The logarithmic slope 
is given by the expression 



a(r) 



-1 



■1 



cs 



1 + f " 1 + cs 

I s 

The mass enclosed within radius r increases as 



M(r) = 47rp s r 



In 1 + 



1 + f 



1 For example. I.Ting and Sut J l)200ft) use the expression p(r 
similar to that advocated by Moore et al. 



(3.6) 



'(r + r s y 



(3.7) 



very 
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or, equivalently, 



M{r) = M 200 g(c) 



ln(l + cs) 



cs 



l + cs 



(3.8) 



Instead of the logarithmic slope (j3.6j) of the radial density profile, in 
order to compare with the results of numerical simulations, we would rather 
use the corresponding value for the mass distribution 



a M (r) = 



d log M(r) 
dlogr 



(3.9) 



We decided to do so because the density profile is prone to significant 
Poisson noise when a narrow radial binning is set, and hence the number of 
particles in each bin is not very high (< 100). Since the numerical evaluation 
of the derivative involved in the computation of a is extremely sensitive to 
these errors, expression (|3.9|) constitutes a much more stable choice to be 
checked against the analytic results. 

In the case of a pure power law, cum = a + 3 = constant. However, the 
NFW profile (jSHJ) yields 



i+— 



In 1 + 



(3.10) 



The limit of this expression as r — > is ajvf(O) = 2, i.e. M(r) oc r 2 . We 
see that, although density diverges as 1/r, the enclosed mass has always a 
finite value, tending to zero when we choose a sphere of null radius. 



Moore et al. 

Based on a set of high resolution simulations, Moore et al. ( 1999bh claim 



that the density profile of dark matter haloes can be fit by the function 

Pm 



p{r) 



(3.11) 



(r/r m )3/2 [l + (r/r m )3/2] 

Th is expression is very similar to the form advocated by .ling and Sutol 
(2000) "or galaxy clusters. The logarithmic slope of this density profile 



a(r) 



1 + 



V r m ) 3 / 2 



reaching an asymptotic value a = — 
The cumulative mass is given by 



l + (r/r m )3/2_ 
1.5 near the centre. 



le is 



(3.12) 



M(r) = M- 



200 



In [1 + (r/r m ) 3 / 2 ] 
In [1 + (r 200 /r m ) 3 / 2 ] 



(3.13) 
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whose logarithmic slope is 



<*M(r) = -- 



(r/r m ) 3 / 2 



(3.14) 



2 [1 + (r/r m ) 3 / 2 ] In [l + (r/r m ) 3 / 2 ] 
implying a growth M(r) oc r 1 ' 5 at small radii. 

3.1.2 Central density 

Obviously, every analytical fit to the mass distribution found in numerical 
simulations must have approximately the same shape. NFW and Moore et 
al. formulae are almost indistinguishable over most part of the halo (r > r s ) 
as long as r m ~ 1.7r s . The only difference between these two profiles is the 
asymptotic behaviour as r — > 0. 

However, we would like to stress at this point that the density profile at 
such small radii is not a prediction of CDM simulations, as it is often stated 
in the literature. The finite number of particles used in numerical N-body 
experiments sets an upper limit to the depth of dark matter potential, and 
the mass distribution at r — > is completely dominated by numerical effects. 
Therefore, any extrapolation of the phenomenological fits to the density 
profile below the resolution limit must be taken with extreme caution, since 
there is no a priori reason why the profile should follow any of the proposed 
formulae at such small distances. 

The minimum radius down to which the mass distribution can be reli- 
ably estimated in numerical simu l ations has been recently investigated by 
Moore et all rtl998al):lKnebe et al.1 (|20001 ): lGhigna et al.M2000h : lKlvDin et"aL 
(|200lh and IPower et alJ ((20021 ) as a function of part icle number, softening 
length or time integration scheme. The first studies by Dubinski and Carlber 
(1991) and NFW contained less than 10 5 particles within the virial radius, 
resolving scales ~ 0.1r2oo- Nowadays, high-resolution numeri cal simulations 
of gala xy clusters attain particle numbers as high as ~ 2 x 10 7 (|Springel et al.l 
2001al ) inside the dark matter halo, which allows to probe the density profile 
at distances smaller than 0.01r2oo- 

(|2002l l. is that the log- 



An important fact, pointed out bv lPower et al 



arithmic slope becomes increasingly shallow inwards, with little sign of ap- 
proaching an asymptotic value at the resolved radii. In that case, the precise 
value of a at some definite cut-off scale (either a fixed physical distance or 
relative to 7"2oo) would not be particu larly meaningful. A similar argument 
has been used by Klvpin et aP (|200lh to explain the different inner slopes 
found by .Ting and Sutol (|2000l ) in terms of the c — M200 relation. 



Observational constraints at galactic scales 

It is somewhat ironic to think that the cold dark matter scenario, originally 
proposed to explain the observed flat rotation curves of spiral galaxies, is 
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nowadays extremely successful in describing the large-scale structure of the 
universe and even the whole process of structure formation, but faces impor- 
tant problems when confronted with the shape of observed rotation curves 
at sub- galactic scales. 

Although rotation curves of nearby galaxies have been routinely mea- 
sured since the early 1970s ( Rubin and Fordl 197fll : Rogstad and ShostakI 
1972h . the analysis of observational data has continued to evolve as larger 
telescopes and improved detect ors became available fo r optical, radio and 
millimetre wavelengths (see e.g. ISofue and R,ubirJ l2001). 

The gravitational mass can be easily inferred from the position-velocity 
diagram, just equating 



GM(r) 



v(r) 



J sys 



sin? 



(3.15) 



where i is the inclination angle and v sys the systemic velocity of the galaxy. 

Combined measurements of HI and Ha or CO emission lines are the best 
tools to probe the dark matter content at galactic scales. Optical rotation 
curves provide high spatial resolution near the centre, while only the neutral 
Hydrogen gas extends far enough in radius to trace the oute r parts. 



O bserved rotat ion curves of dwarf spiral and LSB galaxies dFlores and Primackl 



1994 



Moorelll99llBurkertlll995l:lKravtsov et alJliflQalBorriello and Saluccil 



Salucci 
MM) 



Mill; Ide Blok et ahlEoml: Ide Blok and Bosmalliool iMarchesini et allEoO 
seem to indicate that the shape of the density profile at small scales is sig- 
nificantly shallower than the cusps predicted by both fitting models. 

This discrepancy has been often signalled as a ge nuine crisis of th e 
CDM scenario, and several altern atives, suc h as warm dCohn et alJ feoOO 
Sommer-Larsen and DolgovfeOQll). repulsive dGoodmanl200Ch . fl uid d Peeblei, 
2000h . fuzzy (jHii et all2o7)oTl. decaying dCenl200ll) ann i hilating dKanlinghat et al 
200 0]) or self-inte racting (|Snergel and Steinhardl]l2000l : lYoshida et al.ll2000l : 
Dave et al.ll200ll ) dark matter, have been invoked. 



Unfortunately, it has been proved remarkably hard to esta blish the inner 



slope of the dark mat ter distribution observationally (see e.g.lSwaters et al, 
2002 ) . Some authors dvan den Bosch and Swatersl2001 : Blais-Ouellette et al 
2001 ; I.Timenez et al. 20021 : ISwaters et al.l 2002h claim that a cuspy density 
profile with a < 1 is consistent with current observations, although a shal- 
lower mass distribution is able to explain them as well. Yet, a value as steep 
as a = 1.5 can be confidently ruled out in most cases. 



Observations of galaxy clusters 

Historically, most estimates of the masses of clusters were made from op- 
tical studies of their galaxy dynamics, wherein the motions of individual 
galaxies were used to trace the cluster potentials by virtue of the virial 
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were sen- 



theorem (|Zwickvlll937l : ISmithl l Il936h . Nonetheless, these studies 
sitive to systematic unc ertaint i es du e to velocity anisot r opics, substructure 
and projection effects (jLncevi Il98l iFrenk et al. 1990; Ivan Haarlem et al 
1997). More recent work based on larger galaxy samples and employin; 



199 1). More recent work based on larger galaxy samples and em ploying 
careful selection techniques has lead to significant progress dCarlberget al.l 



carelul selection t echniqu es has lead to sig nifican t prog 
199fil:lden Hartog and Katgertlll996l: iFadda et al.ll l996: 



Mazure et al 



199< 



6; 



TWnni pt. allllMWHCriler nt nl HlfliwHKnnmvi n.nH CrileilbOTlh 

Accurate measures of cluster masses are accomplished via X-ray obser- 
vations of the hot intra-cluster gas and analysis of distortion of background 
galaxies due to gravitational lensing by the cluster mass. 

X-ray mass measurements rely on the assumption that the emitting gas 
which pervades clusters is in hydrostatic equilibrium. Using spherical sym- 
metry and considering purely thermal support, 



1 dP 

p dr 



GM 



(3.16) 



where the total mass profile is determined once t he radial profile s of the 
gas density and temperature are known (see e.g. ISarazinl (jl988bh . The 
gas density can be directly obtained from X-ray images, while temperature 
requires detailed spatially-resolved spectroscopy. The advent of Chandra 
and XMM-Newton satellites allows, for the first time, sufficiently good spa- 
tial and spectral resolution for self-consistent determinations of the density, 
temperature and mass profiles of galaxy clusters. 

In contrast to the aforementioned optical galaxy dispersion and X-ray 
techniques, gravitational lensing offers a method for measuring the projected 
surface density of matter along the line-of-sight which is essentially free from 
assumptions about the co mposition or dynamical state of the gravitating 
material (e.g. Mellier 19991 ). While strong lensing requires compact, dense 
cluster cores and thus probes their central mass distribution, weak lensing 
arclets can be found everywhere across clusters and allow their mass to be 
mapped even at clustercentric distances comparable to the virial radius. 

Clearly, the best approach is to combine the se three methods. The first 
analy sis combining X-ray and lensing studies (jMiralda-Escude and Babull 
1995) suggested that strong lensing masses within 200 kpc of the cluster 
centre might be typically hig her by a fa ctor of ~ 2 — 3 than the X-ray 

e.g.lAlk 



inferred mass. Later work (e.g. lAllenl ll998) highlighted that this discrepancy 
only occurs for non-cooling flow clusters, which typically appear to be more 

dynamically active t han cooling-flow systems. 

B oth X-ray data dEttori et alJ2002bl:IPratt and Arnaudl2002l:IScnmidt et al 



2001 



Allen et all 2001al : lArabadiis et alJl2002l : iLewis et alJ l2002ah . s trong 



(jOguri et al J 12001 ) and weak lensing observations jDahle et al.l l2002h find 
steep slopes consistent with the CD M paradigm (i.e. NFW or Moore et 



al. profiles). However, some authors dTvson et al.lll998l : ISand et alJl2002h 



report evidence of soft cores based on lensing measurements, and some of 
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the X-ray derived profiles are not able to rule out the possibility of a shallow 
(or even null) inner slope. The size of a possible core is constrained to be 
less than ~ 50 kpc, but, for the moment, it seems that current statistics do 
not allow any firm conclusion to be drawn on cluster scales. 



3.1.3 Phase-space structure 



Since CDM obeys the collisionless Boltzmann equation, the velocity disper- 
sion is related to the density field by the Jeans equation. Assuming spherical 
symmetry and neglecti ng bulk velocities (i.e. infall and rotation), the Jeans 
equation reads (see e.g. Binnev and Tremainel 198?! ) 



1 d 

p dr 



{pa 2 r ) + 



GM 



(3.17) 



where a r is the one-dimensional velocity dispersion and (3 is the anisotropy 
parameter, defined as j3 = 1 — cr|/of • If random motions are considered to 
be isotropic {(3 = 0), the radial velocity dispersion is entirely determined by 
the mass distribution. For a NFW profile, 



a 2 Jr) 



1 + 



47rGp s r, 



x 3 (l + x) c " 



ln(l + x) 



l + x 



Ax (3.18) 



from which the corresponding 'universal' profile for a r can be obtained nu- 
merically. 

Following a more empirical approach, Tavlor and Navarro! ( 200 ll ) realised 
that the coarse-grained phase-space density of their galaxy-sized haloes fol- 
lowed a power law 

P 



oc r 



(j 



(3.19) 



with q ~ 1.875 over more than two and a half decades in radius. This 
statement is incompatible with a NFW profile and the Jeans equation, since 
()3.19|) . once substituted into the Jeans equation, fixes both the velocity 
dispersion and the mass distribution of the dark matter halo. The density 
profile thus obtained depends on the normalisation of the velocity dispersion, 
becoming unphysical (decreasing towards the centre) beyond a critical value. 
For this critical value, the density profile is very similar to NFW, although 

the asymptotic slope as r — > is a = 0.7 instead of a = 1. 

Concerning ordered motions, Barnes and Efstathioul ( 198?1 ) found in their 
early work an indication for a 'universal' angular momentum profile, show- 
ing a rough proportionality between j (the specific angular momentum, 
computed in spheri cal shells) and the radial coordinate. More recently, 



Bullock et al.l ( 200ll ) claimed that the cumulative mass with projected an- 



2 When comparing the projected and total angular momentum in cells throughout the 
halo volume, the profiles were found to be very similar (within a factor of less than 2). 
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gular momentum less or equal to j can be described by the function 

M(j) fij 



M 



200 



JO + j 



(3.20) 



where j reaches a maximum value j'200 = Jo/(^ — 1) and £i > 1 is a shape pa- 
rameter. If the angular momentum distribution is assumed to be spherically 
symmetric and monotonically increasing with r, expression (|3.20|) yields 



Jo 



fxM 200 /M{r) 



1 



(3.21) 



Although Bu llock et a l. (2001) claim that the specific angular momen- 
tum distribution in their sample of simulated dark matter haloes is closer 
to being cylindrically symmetric than spherically symmetric, they address 
the question of the behaviour of j(r) averaged over spherical shells, finding 
j oc M or, equivalently, j oc r . 



3.2 Simulation Results 

In order to study the mass distribution in cluster-size haloes, we ran a series 
of high-resolution numerical simulations, using two completely independent 
codes. Since dark matter is the main contributor to the mass budget of 
galaxy clusters, we used the N-body version of the adaptive mesh code 
ART (see Section [2.1. 2 j) . including dark matter only. Furthermore, we also 
simulated the same clusters of galaxies with the Tree-SPH code Gadget 
(Section 12.1. 3JI . which allows us to 

1. Test the consistency of both integration schemes. 

2. Quantify the effects of a baryonic gas on the final density profiles. 

The main properties of our sample of galaxy clusters is briefly sum- 
marised in Section r2.3.1l for both ART and Gadget simulations. Individual 
profiles of all clusters can be found in Appendix lAl 



ART vs Gadget 

Figure l3~T1 shows the ratio between the cumulative dark matter masses found 
in ART and Gadget as a function of clustercentric distance (for details on 
the computation of the centre of mass and radial profiles, see Section |2.2|) . 
Due to the presence of baryons in Gadget, we must rescale the masses by 
the appropriate value of ^cdm = — i n order to compare results: 

AM = M ART / 

m M Gadget \ n m J y ' ' 
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r [h -1 Mpc] 

Figure 3.1: Comparison between ART and Gadget dark matter distributions. The 
relative difference in the cumulative mass is in most haloes (left panel) less 
than 15%. Clusters F, I, J and K (right panel) display bigger discrepancies 
(see text). 

We compute this quantity for each cluster and average the individual 
profiles as discussed in Section 12.2.31 Clusters F, I, Ji and Ki are plotted 
separately on the right panel of Figure 13.11 and were not taken into account 
in the overall comparison. These clusters show discrepancies as high as 40% 
between both codes, but it cannot be attributed to inconsistencies in the 
simulation techniques but to slight differences in timing. Since all these 
clusters are undergoing major merging processes, small differences in the 
time of observation can have dramatic consequences on the morphology 
and dynamical state of the system. Indeed, clusters J2 and K2 are barely 
detectable in ART; in these simulations, we would simply have labelled 
clusters J and K as major mergers (obviously, clusters J 2 and K2 have been 
excluded from these comparison, too). 

As can be seen in the left panel of Figure I3.1| most individual profiles of 
our dark haloes are consistent within ~ 10% accuracy. The mean deviation 
is null for most of the radial range, and it is only slightly biased towards 
larger masses in ART. However, this deviation is small compared to the 
scatter that it is quite likely to be a spurious effect. The fact that ART 
clusters appear to be a little bit heavier both near the centre and beyond 
the virial radius gives further support to this conclusion. 

Baryonic gas 

Another interesting issue is whether baryons affect the shape of the dark 
matter distribution. Non-adiabatic baryonic processes such as cooling, star 
formation and feedback can change the phase-space distribution of the gas 
component, making it gravitationally dominant at the inner parts of galaxy- 
sized haloes. In fact, results from some gasdynamical simulations of galaxy 
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formation seem to indicate that gal actic dark matter density profiles could 
be as steep as p(r) oc r for r — > ( Tissera and Dominguez- Terireirol [l998: 
Thacker and CouchmarJ 200 j ). 



However, most baryons in clusters of galaxies are in the form of a hot 
plasma and their dynamical effects on the ho st dark matter hal o are less im- 
portant even when cooling is considered (e.g. Loken et al.l Eo02). Our results 
seem to support this view, but a word of caution must be said here because 
our simulations take into account adiabatic physics only. In the central re- 
gions of clusters, the typical cooling times are less than a Hubble time, and 
hence a ce rtain amount of cooling is expected to take place. Indeed, cooling 
flows (see Fabian ( 1994^ for a recent review) have been observed in many 
clusters of g alaxies 3 , as well as in numerical simulations including radiative 
cooling (e.g. Pearce et al.ll200oh . 

If baryons played a significant role in the dynamics of the cluster, one 
would expect that they produced an adiabatic contraction of the dark mat- 
ter halo, as it occurs on galactic scales, and hence the mass computed by 
Gadget should be appreciably smaller than that obtained by ART at small 
radii. Since we do not observe such a systematic trend in Figure we con " 
elude that the presence of a hot intracluster gas has no effect on the total 
mass distribution of the CDM component (or, at least, this effect is smaller 
than the 10% scatter we find between the results of simulating the same 
cluster with both codes). 



3.2.1 Universality 

In this section, we address the question of whether the radial dark mat- 
ter density profile of our clusters follows a universal form, comparing the 
sim ulated mass d i stribut ions with the analytical profiles proposed by NFW 
and iMoore et al.l (|1999bh . If the clusters were described by a universal two- 
parameter function, then all the radial profiles should exactly overlap when 
they are properly rescaled. In the case of NFW formula (j3.2j) . the scale 
factors are the characteristic density p s in one axis and the characteristic 
radius r s in the other. 



Fitting procedure 

The values of p s and r s are obtained by a x 2 fit to the cumulative dark 
matter mass. As we discussed in Section [3.1.11 the density profile is much 
noisier than the cumulative mass, and hence much more sensitive to the 
chosen radial binning. We took equally-spaced logarithmic bins in order 

Nevertheless, high-resolution Chandra images and XMM-Newton spectra sh ow that 
the p hysics of the ICM within the central 100 kpc can be extremely complex jFabianl 
I2002TI , and it is likely that some form of heating reduces the mass cooling rates by a factor 
up to ten. 
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Table 3.1: Best-fit values of NFW parameters p s (10 15 /i 2 M Mpc~ 3 ), r s (h^ 1 kpc) 
and c = r2oo/fs- Merging haloes are typically more extended than relaxed 
ones, and they are less accurately fit by a NFW profile. 



to give more statistical weight to the inner regions, considering only radii 
between 0.05 r2oo an d r^uo ^° cnm pute the r.m.s deviations from the fit. 
In the less massive haloes, our resolution limit of 200 dark matter particles 
(plus 100 gas particles in the case of Gadget simulations) placed a stronger 
constraint on the lower limit. In total, 26 points have been used for the fit 
in all clusters. 

Actually, we fitted the mass distribution to expression (|3.8|) in terms of 
r2oo and c, and then p s and r s were computed from the best values of these 
parameters. We imposed the physical constraints that c has been let to vary 
from c = 1 to c = 15, while a first guess for r2oo was obtained from the mass 
profile itself (see Table l2~2*)) . The best-fit value of this quantity was allowed 
to fluctuate between 0.8 and 1.2 times the initial guess. 

Table l3~Tl summaries the values obtained by this method. The accuracy 
of each profile is expressed as the root mean square deviation of the numerical 
profiles from the analytical fit. For relaxed clusters, this quantity is always 
within a few percent, the quality of the fit being slightly better for more 
massive clusters. Minor mergers deviate about 6% and the discrepancy rises 
up to more than 10% in the major mergers. 

At first sight, the well known c— M200 relation fe.g. lNavarro et al.lll997l ) 
is not apparent in this table. The scatter in the concentration parameter 




Figure 3.2: Top panel: dark matter density profiles (squares) scaled by the charac- 
teristic density and radius of each halo. Red lines indicate NFW fit; blue 
is used for Moore et al. Bottom panel: logarithmic slope au of the mass 
profile. Merging clusters display a pure power-law behaviour near the 
centre, while the slope of relaxed structures decreases gently with radius. 



seem s to be quite large in the mass range covered by our sample (see also lJind 
2000). However, there is a trend in the sense that c seems to be correlated 
with the dynamical state of the cluster. A more detailed discussion on 
this issue is given below, in connection with the effects of merging and 
environment. 



Mass profiles 

We plot the radial density profiles of our clusters in the top panel of Fig- 
ure rescaled by the best-fit values of p s and r s given in Table IXT1 The 
phenomenological formulae proposed by NFW (red) and Moore et al. (blue) 
are also shown for comparison. Since we have fitted our dark matter haloes 
to the NFW form, we have assumed r m = 1.7r s and p m = 6p s to represent 
the Moore et al. 'universal' profile in the figure. 

We have plotted ART/Gadget data as empty/solid squares, respec- 
tively, and the innermost bin is set by the constraints explained in Sec- 
tion 12.2.11 The formal resolution is similar in both codes: the finest grid 
level in ART and the gravitational softening length in Gadget were set to 
2 h- 1 kpc. 

From the top panel of Figure I3.2| we learn that most of the clusters of 
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galaxies in the present sample can be reasonably described by either NFW 
or Moore et al. analytical fits, although some of them (specially major 
mergers) display strong deviations from the 'universal' profile for r > r2oo 
(a few times r s ; see table 13.1]) . 

Nevertheless, a log-log plot of the density profile offers very limited infor- 
mation about the subtleties of the radial distribution (in fact, any deviation 
smaller than a factor of 2 is not easily perceived, given the scale of the fig- 
ure). In order to study the shape of the mass distribution in the innermost 
regions, we plot the local value of the logarithmic slope of the cumulative 
mass profile au = dlogM/dlog function of radius (bottom panel of 

Figure |2HJ)- The difference in the logarithmic slopes predicted by NFW and 
Moore et al. is mostly evident for r < r s . 

We see that relaxed and merging clusters show a remarkably different 
behaviour near the origin. Focusing only on the relaxed subset of the sam- 
ple, we find that NFW provides a good fit to our data. The steeper slope 
predicted by Moore et al. formula can be confidently ruled out at radii as 
large as one order of magnitude above our res olution limit. 

However, in agreement with the results of Power et al. ( 20021 s ). we find 



that ttj,/ does not tend to any asymptotic value at all, but steadily in- 
creases towards the centre. Even more resolution is required in order to 
reach smaller radii (at the kpc scale) and test whether the logarithmic slope 
stabilises at a constant value of 2 (as in NFW fit) or, on the contrary, it still 
increases until Mar 3 (i.e. a constant density core). 

This seems indeed to be the case when the simulated profiles are plotted 
slightly beyond the resolution limit. Yet, we would like to stress that the 
increasing slope at r < 0.1r s may well be a numerical artifact induced by 
poor resolution. A this point, we prefer to be conservative and wait for 
higher-resolution data before reaching any conclusion. 



Dependence on environment 

Relaxed clusters seem to follow indeed a universal profile dependent on two 
scale parameters only (albeit a small scatter exists for r > r2oo)- This profile 
can be described with reasonable accuracy by NFW expression up to our 
resolution limit, being the innermost values of the logarithmic slope cum 
clearly inconsistent with the proposal of Moore et al. 

However, the situation is quite different for merging clusters. Minor 
mergers feature logarithmic slopes systematically steeper than those of re- 
laxed systems. The shape of their radial density profiles is not exactly NFW 
or Moore et al.-like, but lies somewhere in between both fits, in many occa- 
sions closer to the latter (a« = 1.5) and sometimes even steeper. 

Major mergers are an extreme case. Although they can be roughly ap- 
proximated by NFW or Moore et al. formulae, significant deviations occur 
throughout the whole cluster. Nevertheless, our results hint that these re- 
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Figure 3.3: Relation between M200 and concentration parameter c (left panel) and 
characteristic density and radius (right panel). Black squares represent 
relaxed systems, empty squares show minor merg ers a n d star s are used for 
major mergers. Observational data from ISato et a 1. 1 l|200Ctl is indicated 
by the dashed line. 



cent mergers share some similar characteristics, the most remarkable being 
a power-law profile for more than one decade in radius. At the inner part, 
the logarithmic slope remains approximately constant. The exact value of 
olm shows some scatter, but the average seems to be slightly lower than 1.5. 

A conjecture that deserves further investigation is that violent relaxation 
and subsequent equipartition of energy drive the inner part of the system 
into an isothermal state, and thus olm = 1 immediately after the merger. 
The cluster would then slowly relax, and the density profile would become 
increasingly shallow with time until virial equ ilibrium is reached. This wo uld 
be consistent with the scenario proposed by S alvador-Sole et al. I l|l99rf ). m 



which the halo formation time is defined by the last major merger. 

The fact that major mergers rearrange the inner structure of dark matter 
haloes should leave some common imprint on the radial mass distribution. 
Apart from the steeper slope of the density profile, one of the features shared 
by all the merging systems in our sample (already shown in Table l3~T|) was 
that merging haloes tend to be much less concent rated than relaxed ones. 
This is in agreement with the results of Jin j ( 2f)0Cl ) , who concluded that the 
less virialised haloes (but still fitted by a NFW profile, similar to our 'minor 
mergers'), have a mean concentration about 15% smaller than systems in 
virial equilibrium. 

In Figure ESI the c — M200 (or, equivalently, r s — p s ) relation is shown 
for our cluster sample. We see a clear dichotomy between the behaviour 
of the c — M200 relation of relaxed and merging systems. While relaxed 
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Figure 3.4: Accuracy of NFW fit for the simulated clusters 



clusters follow the usual trend of bigger concentration for smaller mass, 
major mergers are so extended that their concentrations are much lower 
than those of the most massive systems. 

Due to the large variations that dark matter haloes experience during a 
merging event, it seems unlikely that a c — M200 relation can be established 
during these transient stages of their evolution. A similar situation arises 
in minor mergers, where the effects on the concentration (on the whole 
density profile, indeed) depend strongly on the specifics of the encounter, 
such as the mass ratio of the CDM haloes, their relative velocity, their former 
concentrations, etc. 

The relation seems to be a little tighter in the r s — p a plane, shown on 
the right panel of Figure EZD Although major mergers still show a different 
behaviour than relaxed clusters, their profiles have both larger values of the 
characteristic radius and lower central densities. This trend agrees well with 
recent X-ray obser yational estimate s of the mass profiles of objects between 
10 12 and 10 15 M (jSato et al .11200(1 dashed line in Figure ESI). 



Accuracy of NFW fit 

One way to quantify to what extent do clusters follow a given 'universal' 
profile is plotting the deviations with respect to the proposed form. Since we 
have fitted our CDM haloes by a NFW function in order to scale both axis 
in Figure E21 we will present our results taking this profile as a reference. 

We saw in Table 13.11 that the r.m.s values obtained for the NFW fit 
indicated accuracies of the order of 5 — 10%. However, we recall that this 
only applies for the cluster region that has been fitted (i.e. from 0.05 r2oo 
to r2oo)- Outside this area, the profile of our dark matter haloes is not so 
well described by the NFW form. 

The r.m.s deviation for the whole fit offers only a limited test of the 
ability of NFW profile to trace the mass distribution found in numerical 
experiments. We can obtain more information from Figure 13. 4( where the 
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relative difference between this profile and the cumulative mass of our simu- 
lated haloes is plotted against clustercentric radius. Excluding major merg- 
ers, our clusters do not deviate appreciably from the NFW form throughout 
most of the radial range (0.1 — l)r2oo- Even in major mergers, the deviation 
is seldom larger than 20% in this region, although the formula proposed by 
NFW offers a poor description of the shape (i.e. the logarithmic slope) of 

the mass distribution in these systems. 

However, as pointed out by many authors (e.g. Moore et al.lll9 99b. 



sec 



Section I3.1.1)) . departures from the NFW formula seem to be similar for 
all clusters, displaying a characteristic 'S' shape. Nonetheless, the mass 
excess detected at very small radii does not necessarily imply a steeper 
inner slope of the radial de nsity profile. For example, the form proposed by 



Tavlor and Navarrol ()200ll ) (see Section I3.1.3|) results in more concentrated 
CDM distributions, but the asymptotic slope a(r —* 0) ~ 0.7. Yet, the 
fact that our simulated dark matter haloes show systematic deviations with 
respect to the NFW formula hints that it is not a good approximation to 
the shape of their density profile at the innermost regions. 



3.2.2 Phase-space structure 

In addition to the radial mass distribution, the kinetic structure of simulated 
CDM haloes can offer interesting insights into the formation of galaxies and 
galaxy clusters. On cluster scales, rotation is expected to give a negligible 
amount of support against gravitational collapse because the angular mo- 
mentum of the system is dominated by the contribution of random motions. 

In this section, we will discuss the phase-density structure of our clusters, 
separating the global and random compo nents and comparing our result s 
with previous work bv lBullock et al. ( 200 il l and Tavlor and Navarrol ( 2001 ). 



respectively (see Section f3.1.3Jl . The details concerning the computation 
of the angular momentum and velocity dispersion profiles can be found in 
Section 



Velocity Dispersion 

As was discussed in Section 13.1. 3| the velocity dispersion profile of a collis- 
sionless dark matter halo is related to the total mass distribution by virtue 
of the Jeans equation (|3.17|) . The average radial profile of the spherically- 
averaged phase-space density p/a 3 has been computed for the whole sample 
of galaxy clusters, and it is plotted in Figure 1531 Individual profiles display 
a remarkably low scatter, particularly if we take into account that 

1. Major mergers have been included in the average. 

2. The profiles have not been normalised. 
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Figure 3.5: Phase space density profiles for our h aloes (dots) compared to the fit 
proposed bv lTavlor and Navarro! (|200lft (p/a 3 oc r -1 - 879 ) 



The phase-space density of our simulated clusters is very well fitted by 
the simple power-law 



p/a 3 



2.24 x 10 4 



-1.875 



(3.23) 



where p, a and r are expressed in units of h 2 M Mpc~ 3 , km s _1 and hr x 
Mpc, respectively. The slope of the power l aw (-1 .875) corresponds to the 
self-similar solution derived bv iBertschinger (1985) for secondary infall onto 
a spherical perturbation in an otherwise homogeneous Einstein-de Sitter 
universe, and it has been found to describe also extre mely well the phase- 
space density structure of galactic-size CDM haloes ( Tavlor and Navarro! 

An interesting issue which would be worth to investigate is whether 
the normalisation depends somehow on the scale of the halo, the assumed 
cosmological model, or any other factor. Un fortunately the phase-space 
profiles computed by Tavlor and Navarro! ( 200lh have been scaled arbitrarily, 
so it is impossible to compare with our result (|3.23j) . 

This correlation between density and velocity dispersion profiles has been 
used by Tavlor and Navarro! (j200ll ) to infe r the radial mas s distribution of 
the dark matter from the Jeans equation. Hiotelisl ( 2002bh has generalised 
this solution for the case of anisotropic velocity tensors, using an approxi- 
mation for the anisotropy parameter 



P(r) =Pi + 2f3 2 



r/r* 



1 + {r/r^f 



(3.24) 
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Figure 3.6: Anisotropy parameter (3. See text for description of curves 



where B\ and 8?. are free parameters, and r* is twice the virial radius 
( Carlberg et alilS IColfn et aljEonfl) . 

In Figure 13.61 the anisotropy profile of our simulated clus ters is com- 
pared w ith expression (|3.24|) . using the best-fit values found bv lCohn et al 
(20003), 8\ = 0.15, and 82 = 0.5 (dashed lines). Although the average pro- 



file of the anisotropy parameter seems to follow a similar trend in all haloes, 
the assumption of a 'universal' anisotropy profile seems hardly justifiable 
given the scatter of each individual profile around the mean value 4 . 

In any case, the simple formula (|3.24|) proposed by Carlberg et al. ( 1997h 
does not provide a good description of the shape of the anisotropy profiles of 
our cluster sample. The exception could be the clusters with major mergers, 
which are much closer to isotropy in velocities than the rest, but we would 
need to modify the values of 8\ and 82 and change r* by r2oo- Nevertheless, 
expression (|3.24[) is a much better approximation to our results than purely 
isotropic profiles {8 = 0, dotted line), which can be ruled out in relaxed 
clusters and minor mergers at a several sigma confidence level over a broad 
radial range. 

Our numerical results can be better fit by the function (solid lines in 
Figure EEJ) 

r/r* 



8{r) = & 



1 + r/r* 



(3.25) 



where /3 max = (0.7,0.6,0.5) and r* = (0.07, 0.1, 0.4)r2oo for relaxed clusters, 
minor and major mergers respectively. We insist, nonetheless, that this fit 
is a poor approximation to the radial dependence of 8 in any individual 
halo, and it could only be applied in a statistical se nse, such as in the 
determination of the average mass profile attempted hv lHiotelisI (j2002hl V 



4 We note, however, that the numerical derivation of j3 is not straightforward, since this 
parameter is extremely sensitive to inaccuracies in the velocity profiles, particularly near 
the centre. 
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Figure 3.7: Angular momentum in the simulated haloes. Left panel: Contribution of 
bulk rotation (solid squares) and random motions (empty squares), jar 
is drawn as a dashed line. Right panel: Eccentricity of the individual 
orbits (haloes I and J2 are shown separately). 



Angular momentum 

The origin of angular momentum in dark matter haloes can be understood in 
terms of lineal tidal torque theory, in which prot ohaloes acquire their angular 
momentum because of the surrounding shear (|Peebles!ll96fll : Doroshkevichl 
197(1 IWhitelll984l : IPorciani et aLll2002al lbh . 



In the left panel of Figure 13.71 we plot the angular momentum of our 
haloes computed in spherical shells around the centre of mass. Solid squares 
indicate the mean value of j in each shell (i.e. rigid body rotation). The 
amount of angular momentum due to random motions (i.e. tangential ve- 
locity dispersion) is represented by the open symbols. 

The first thing that can be appreciated in Figure 13.71 is that velocity 
dispersion is the dominant source of kinetic energy in the system. The 
amount of angular momentum contributed by random motions is higher 
by approximately one order of magnitude than that corresponding to bulk 
rotation. 

The relative scatter in the radial profiles (when normalised on both axes 
as indicated in the figure) is larger for the mean angular momentum than 
for the r.m.s. value. Near the centre, the situation is less clear, since at 
very small radii the low number of particles (and possible misidentifications 
of the centre of mass and velocities) makes difficult to disentangle ordered 
from random motions, as well as radial and tangential direction s. 

Building upon the pioneering work of lBarnes and Fffitathloiil l|l987h . the 



question of a universal angul ar momentum profile has been recently investi- 
gated by several authors (e.g. Bullock et al.ll20fll ; van den Bosch et al.ll2002l : 
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Chen and ,Tin"gll2002h . The spin profiles found in n umerical simulation s grow 



approximately linearly with radius (or as j oc r 11 . lBullock et aL 2001 '). Our 
results are consistent with this behaviour (dashed line in the left panel of 
Figure ETT)) . but they deviate from this simple relation near the centre. How- 
ever, the discrepancy (~ 2 over j oc r at r ~ 0.01r2oo) is comparable to the 
scatter in the profiles. We must take into account that, at these radii, the 
lower number of particles may well lead to an overestimate of the mean 
angular momentum with respect to the random component. 

Once the angular momentum distribution of the dark matter is known, 
it is interesting to study the structure of the orbits of individual particles 
around the centre of mass. More precisely, the eccentricity of these orbits has 
important dynamical consequences, since it controls directly the maximum 
depth that a dark matter particle is able to reach within the gravitational 
potential of the dark matter halo. The ideal case e = corresponds to a 
configuration in which all particles would orbit at constant radius with no 
mixing between concentric shells, while e = 1 (often assumed in spherical 
infall models) implies that every dark matter particle goes through the very 
centre of mass of the cluster. In general, large values of the eccentricity lead 
to a situation in which the core is mainly composed of particles coming from 
the outer shells of the halo that happen to be near the pericenter of their 
orbits. 

We plot the the eccentricity profiles of our haloes as a function of the 
apocentric radius r m in the right panel of Figure 13.71 Both quantities have 



been computed according to the prescriptions given in Section [2231 Haloes 
I and J2 are plotted separately as dashed line, since they feature much lower 
eccentricities than the rest of the sample. Except for these clusters, the 
rest of the sample exhibits a similar profile (as can be seen from the small 
scatter). For large apocentres, the eccentricities of different clusters show 
different behaviours, although the trend is that orbits are in general more 
radial as we move out to the turn-around radius. 

Despite the small dependence of e on r seen in Figure l3~Tl we will make 
the approximation e ~ 0.5 (indicated by the horizontal dotted line) in order 
to model the orbits of dark matter particles. With this value, the minimum 
radius reached in the orbit is one third of the apocentric distance, while it 
would be ~ 0.43 r m for e = 0.4 and 0.25 r m for e = 0.6. Since the error we 
make is not large, we will use this simple approximation in our treatment of 
angular momentum in the spherical collapse formalism (Section I3.3.3|) . 



3.3 On the physical origin of dark matter profiles 

From the theoretical point of view, a number of plausible arguments have 
been advanced to try to explain the dynamical structure of dark matter 
haloes. The controversy regarding the 'universal' density profile and the 
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asymptotic value of its logarithmic slope at the centre has stimulated a 
great deal of analytical work attempting to find the physical reasons behind 
this similarity in the mass distribution, as well as possible correlations with 
halo size, environment, underlying power spectrum or even the nature of 
dark matter particles. 

The structure of collapsed and virialised objects such as galaxies and 
clusters poses a real challenge to our understanding of structure formation 
in the universe. The basic problem of the collisionless collapse of a spher- 
ical perturbation in a n expanding backgroun d w as addressed first by the 
two seminal papers of lGunn and Gott jl972h and lGunnl jl977h . where the 
cosmological expansion and the role of adiabatic invariance were first intro- 
duced in the context of the formation o f individual objects. 

The next step w as accomplished by Fillmore and Goldreichl (jl98 4) and 
Bert schinger (13, who found analytical predictions for the density of col- 



lapsed objects seeded by scal e - free primordial perturbations in a flat uni 



verse. 



Hoffman and ShahamI ()l985l ) generalised these solutions to realis- 



tic initial conditions in flat as well as open Friedmann models. Modifica- 
tions of the self-similar collapse model to include more realistic dynamics 
of the growth process have b een attempt e d (e.g. lAvila-Reese et a~ _ 19981: 



Henriksen and Widrowl ll 9991: l^okad 1200(1 IKuIII Il999l: ISuhra.manian et al l 
2000). Several authors (e.g. ISver and Whitelll998| : ISalvadoi^ Sole et al.lll3 



Nusser and Shethlll999l ) argue that the central density profile is linked to the 



accretion and merging history of dark matter substructure. 

In this section, we will focus on the predictions of spherical collapse 
theory, showing that this simple model is able to explain the general fea- 
tures of the mass distribution of our simulated haloes described in previous 
sections. First we explain our numerical implementation of the spherical 
collapse model and the primordial initial conditions. We then compare the 
resulting analytical profiles with our simulations. Finally, a brief discussion 
on the role of hierarchical merging is also included. 



3.3.1 Spherical collapse 

The most simple way of addressing the problem of structure formation is 
to assume spherical symmetry. For a homogeneous and isotropic universe, 
General Relativity leads to the Friedmann equation 1)1. 10|) . Since we are in- 
terested on scales much smaller than the horizon, we can restrict our treat- 
ment to Newtonian dynamics. In this case, the equation of motion for a 
Lagrangian shell of mass M can be derived from energy conservation 

t(r) _£W_g_ G ^(r) + ^_ (r) ( 

m 2 r 

where M(r) = Mj(r,) is the mass contained within the shell, p\ is the 
vacuum energy density, and the subscript 'i' denotes initial conditions. If 
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we assume the universe to be homogeneous, 



4-7T ■ ■ q 
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Substituting in (pO^l . 



r 2 ^4vr , f&r? + A r 3 _ (H iU ) 2 4vr . H^rf + 0\r? 



3 F 



3ff 2 



Since = we can make a coordinate transformation 

r(t) = na(t) 



(3.27) 



(3.28) 



(3.29) 



and recover the Friedmann equation simply by multiplying both sides of 
(EE2ED by a factor : 



a 

Hi 



ft^a' 1 - tt\a 2 = 1 - fi^ - n\ 



Now, we introduce a small mass perturbation 
a / n \ / \ M(n) 

ly % > 1 y 11 47TO* rf-T^ 

3 i 'mPc' i 



1 



(3.30) 



(3.31) 



at an early epoch m, so that the growth is still in the linear phase. The 
initial positions and velocities of the spherical shells of matter are slightly 
modified, keeping only terms to first order in Ajfa): 

ri^rifl-^Aiin)) ; ^ ~ %, I 1 - ?A,(r.) j (3.32) 



and the initial specific energy of each shell is thus 



(H iri ) 2 
2 

(H iri ) 2 



2A, 



A, 



n\ I i 



A; 



A, 



(l-0^-OX)--^(4 + 0^-20X) 



Since a; >C 1, 



f2 m a, + Qa + O k a, 



~ 1 ; ft} 



n f 



f2 m a,- + + f2 k a. 



and energy conservation (J3.26|) leads to the equation of motion 



2e ?: 5 

^ --A,- 



(tf^) 2 - 3 1 H 2 
valid for a spherical shell enclosing a mass Mj. 



^a" 1 - n\a 2 



(3.33) 

~ 
(3.34) 

(3.35) 
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Before turn-around 

According to (|3.35j) . the evolution of a single spherical shell would thus be 
similar to that of a Friedmann universe. For high enough values of A$, 
the shell reaches a maximum radius r m at a turn- around time t m and re- 
collapses again to a point singularity. For an Einstein-de Sitter universe, 
equation (j3.35|) can be solved analytically, giving 



This is also a valid approximation for shells reaching the turn-around 
radius r m before the cosmological constant term starts to dominate the 
expansion. Since the shells need at least another t m to virialise, in the 
ACDM model discussed in Section ITTT1 expression (|3.36[) can be applied to 
estimate the maximum expansion radius and time for the inner part of a 
virialised cluster. For the outer shells, however, the equation of motion 
()3.35j) must be integrated numerically in order to find the trajectory up to 
the maximum radius. 

In absence of shell-crossing, gravity would make the shell collapse back in 
a time-symmetric motion, so that the shell will reach the origin at T = 2t m . 
Since shells are assumed to be composed of collisionless cold dark matter 
particles, they will simply pass through the centre, describing an oscillatory 
motion with amplitude r m and period T. 

After turn-around 

Equation (|3.35|) holds as long as the enclosed mass remains constant. How- 
ever, as a shell re-collapses, it will cross the orbits of the inner shells that 
are bouncing around the centre. When this happens, the assumption of 
constant mass breaks down because 

1. The mass enclosed by the outer shell decreases 

2. The mass enclosed by the inner shells increases 

To take into account the first point above, our model assumes that the 
final mass distribution can be approximated locally by a pure power law. 
Since dark matter particles are expected to spend most of the time in the 
outermost regions of their orbits (see also Section 13.3.31 on angular momen- 
tum), we only need this approximation to be valid over a relatively narrow 
range below the maximum radius. Near r m , the enclosed mass is taken to 
vary with r as 



3n 



7T 



(3.36) 



r. 



5Ai 



2^(5^/3)3/2 




(3.37) 



where aj,;(r m ) is the local value of the logarithmic slope. 
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At first sight, this might seem similar to the classical approach followed i n 
early works based on self-similar spherical collapse (e.g. Bertschinger 1985), 
but in those cases the dark matter distribution was indeed assumed to be 
a power law, whereas in our model this ansatz is only an approximation to 
compute the orbits of particles in a given shell. The final mass profile is 
computed self-consistently as a sum of contributions from all shells, which 
have different values of ajw(r m ). 

After turn-around, the particles belonging to a shell will oscillate (or, 
more generally, orbit) in the gravitational potential created by the dark 
matter halo. Taking the ansatz ()3.37j) . the gravitational potential is approx- 
imated as 



$(r) - $(r m ) 



GM(r) 



GMi 



dr 



r^ocMirm) - 1] 



r 



(3.38) 



The probability of finding a particle inside radius r is just proportional 
to the fraction of time that it spends within r: 



1 r dx If 
P(r,r m ) = - -rr = r / 
hn JO u r yx) i m Jq 



dx 



v r (x) t m J y/$(r m ) - 



(3.39) 



Note that, in this form, expression ()3.39j) is valid regardless of the as- 
sumed ansatz for M(r). Taking different prescriptions to compute the po- 
tential, such as an isothermal profile (ajvf( r m) = 1 f° r every shell, so <3?(r) is 
actually logarithmic) or considering M constant (ajvf = 0, Keplerian poten- 
tial $(r) = —GM/r) does not lead to significant variations in the probability 
P(r,r m ), and hence in the resulting mass distribution. 

If phase mixing is considered to be efficient, particles initially on the same 
shell will be spread out from r = to r = r m a short time after t m . For 
the sake of simplicity, we will consider that phase mixing is instantaneous, 
so immediately after turn-around the shell is transformed into a density 
distribution whose cumulative mass is proportional to P(r,r m ). 

Effect (2) listed above has dramatic influence on the final structure of the 
dark matter halo. If the recently accreted shells have a non-zero contribution 
to the mass within r (i.e. P(r,r m )dM), the mass inside the inner shells 
whose maximum radius was r changes from Mj(r) to 



M(r) = Mi(r) + M add (r) = K(r)Mi(r) 



(3.40) 



where M ai w( r) accounts for particles belo nging to outer shells. To compute 
M adc i( r ) (see lZaroubi and Hoffman] 19931 ). we must integrate the contribu- 
tion of every shell whose maximum radius is larger than r, up to the current 
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turn-around radius R m : 



M add (r) = J —±± P{r,x) dx (3.41) 



Gunnl (Il977h was the first to apply the concept of adiabatic 



mvariance 



to the secondary infall problem. If the orbital period of the inner particles 
is much shorter than the collapse timescale of the outer shells, the dynamics 
of the inner shell admits an adiabatic invariant 

J r = -!- <j> v r (x) dx (3.42) 

also known as the radial action. For a potential of the form (|3.38j) . the radial 
action J r is proportional to W r m M(r m ). Therefore, the maximum radius 
r m of the inner shell decreases to a final value given by the implicit equation 
r = F(r)r m , where 

F( ) = = Ml = Mi (3 43) 

U K(r) Mi + Afadd(r) Mi + Mad d (F(r)r m ) { ' > 

and whose solution must be obtained numerically for each shell. 
Summary 

To summarise, the numerical procedure to compute the final radius of a 
Lagrangian shell of matter involves the following steps: 

1. Set the initial mass Mj according to (|3,27|) 

2. Integrate the equation of motion (|3.35|) up to the maximum radius v m 

3. If t m > to 

• The shell is still expanding, r = r(to) 

4. If t m < t 

• Solve (|3.43|) to compute the 

• Add the contribution P(r',r)dM to M add (r') 

5. Repeat the process for the next shell (closer to the centre). 
3.3.2 Initial conditions 

The spherical collapse model allows us to compute the mass distribution 
arising from a primordial fluctuation Aj(rj), but does not say anything about 
the shape of this function or its physical origin. Nevertheless, it is important 
to note that the final density profile is entirely determined by this initial 
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condition. In the spherical collapse paradigm, the case for universality in 
the radial profiles of dark matter haloes can be reformulated in terms of 
uni versality in the primordial de nsity fluctuations that fix Aj(rj). 

Hoff man and Shahaml (|l985l ) suggested that, according to the hierar- 
chical scenario of structure formation, haloes should collapse around local 
maxima of the smoothed density field. The statistics of peaks in a Gaus- 
sian random f i eld ha s been extensively studied in the classical paper by 



Bardeen et al.l (|l986i .hereafter BBKS). One of the best known results of 
BBKS is the expression for the radial density profile of a fluctuation centred 
on a primordial peak of arbitrary height: 



CO 



7(1 -7 2 ) 



R 2 

l 2 ^(r) + -^V 2 V>(r) 



(3.44) 



where do = C(O) 1 / 2 is the rms density fluctuation, ip(r) = £( r )/°"o is the 
normalised two-point correlation function, vctq is the height of the peak, 
and 7 = o" 2 /(<T2<7o) and i?* = V^oi/o^ are related to the moments of the 
power spectrum 

/ P(jfe) fc 2 ^' +1 ) dk (3.45) 
to 

The function #(7, 71/) parameterises the second derivative of the density 
fluctuation near the peak. BBKS suggest the approximate formula 
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(3.46) 

valid to 1% accuracy in the range 0.4 < 7 < 0.7 and 1 < 7^ < 3, which is 

the scale relevant for both galaxies and galaxy clusters. 

Expression (13.441) is o ften quoted in the literature 5 (e.g. lHoffmanl ll988: 



Lokas and Hofrman^OOnl : IPel PopoIo et alJEnionl : Iffiote lis 2002a|) as the ini- 
tial condition Aj(rj), even though BBKS explicitly warn in their paper 

"Unless the filter is physical (as in pancake models), these profiles 
cannot be used for hydrodynamic or N-body studies of collapse 
since substructure must be included" . 

Although BBKS referred to a different problem (see our discussion in 
Section 13.3.41 regarding the effects of merging and substructure) , the point 
we want to stress is that (|3.44|l describes the smoothed density profile around 
a local maximum of the smoothed density field, 



6 f (r)= W f {r-x)5{x)d 6 x 



(3.47) 



5 Some authors (Hoffman, private communication) use indeed the prescription described 
below. Yet 1)3.44^ is stated in their paper to determine the initial conditions. 
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where the function Wf(f—x) is a smoothing kernel that depends on a certain 
filtering scale Rf. 

However, the smoothed density profile ((5/(r)), given by equation (|3.44j) . 
is never equal to the mean value of the local overdensity (<5(r)} that we need 
to integrate in order to compute Aj(rj): 

Ai(n)=4TT {5(x))x 2 dx (3.48) 
J o 

Comparing this expression with (|3.47j) . we see that A^(r^) is equivalent 
to (5/(0) as long as Wf is taken to be a spherical top hat function of radius r,. 
In phase space, the Fourier transform of this function results in significant 
oscillations near the cut-off scale, so it is better to choose a Gaussian filter 

with Rf = 0.64rj, such that the enclosed mass is approximately the same as 
in (ETiHl) . 

We are interested in the physical density profile around a local maximum 
of the smoothed density field. First, we fix the scale of the fluctuation Rf. 
Then, according to (|3.48j) . Aj(rj) = (5 r (0), where (5 r (0) is the central density 
smoothed on a scale ro = 0.64rj. Furthermore, we want to impose the 
condition that r = is a maximum of the density field <5/(r) smoothed on a 
scale Rf. BBKS show that the probability distribution of (5 r (0) is a Gaussian 
with mean 



for ro > Rf and 



<dr(o)> = vf-^ - li^m^i^ i _ (3.50) 

J o-of 1 - 7/ <TQf \ <ri h off J 



(5 (0)) = Vf ?k (3.51) 

O-Qf 



for ro < Rf. The moments 

°i* = 2^/ P*(k)k^ +1 Uk (3.52) 

are computed from 

Pf(k) = PACDM(k) exp [-(kRf) 2 } (3.53) 

and 

P h (k) = P\cDM(k) exp 



-k 2 



m + r 



(3.54) 



Thus, instead of (|3.44j) . we have used ()3.50|) and ()3.51|) to set the initial 
condition Aj(rj), assuming that the cluster collapses from a Vf&Qf peak on 
scale Rf. 
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3.3.3 Angular momentum 

Two d ecades after the seminal paper by Gunn and GottI ( 19721 ) , White and Zaritskvl 
dl992l 'l introduced the idea (although see 6 also the pioneering work of Gurevich and Zvbinl 
1988) that angular momentum prevents dark matter particles from reaching 
the origin. 

For pure radial orbits, the mass in the centre is dominated by the con- 
tribution Madd from the outer shells when the profile at turn-around is 
shallower than isothermal ( Fillmore and Goldreichl 19841 '! . The final density 
distribution found by these authors w as p(r) oc r~ 2 , ind ependent on the ini- 
tial lo garithmic slope. More recently, Lokad ( 2000) and lLokas and Hoffman! 
(2000) found a similar result consid ering non-self-sim i lar p rimordial fluctua- 
tions based on the peak formalism. Del Popol o et al. (2000) obtains slightly 
shallower profiles, but they rely on IWhite and Zaritskvl jlflflj ) to override 
the 'crowding effect' associated to radial orbits. 

Angular momentum arises in linear theory ()White!ll984h from the mis- 
alignment between the asymmetry of the collapsing region (i.e. the inertia 
tensor) and the t idal forces it exper iences during the lineal regime. Also, vi- 
olent relaxation (|Lvnden- Bell 1967) transforms the radial infall energy into 
random velocity dispersion, due to two-body encounters taking place during 
shell crossing. 

The total amount of angular momentum acquired by the system is, how- 
ever, an open questio n. The usual approach (jAvila-Beese et al.lll998l : lNussei 
2001 ; Hiotelid 2002al ) consists in assigning a specific angular momentum at 
turn-around 

(3.55) 



J oc 



^GMr 



With this prescription, the orbit eccentricity e is the same for all particles 
in the halo ( Nusserl 1200 a \ . As can be seen in Figure 13.71 this is a fair 



approximation for all but two of our simulated clusters. Therefore, we will 
assume hereafter a constant eccentricity e = 0.5 for the orbits of every dark 
matter particle after turn-around. For this value of the eccentricity, the 
pericenter radius of the orbit will be given by 



1 



1 + e 



e rn 
''max — 



(3.56) 



where r max is computed from adiabatic invariance ()3.43[) . following the pro- 
cedure explained in Section 13.3.11 In this case, angular momentum adds 
the usual effective term j 2 /(2r 2 ) to the gravitational potential 3>(r) to be 
substituted in equation (|3.39|) . Although this changes the actual value of 
the radial action, J r is still proportional to y r m M(r m ). Therefore ()3.43|) 
still can be used to compute the collapse factor F(r). The assumption of 
spherical symmetry leads to angular momentum conservation, which implies 
constant e during the contraction. 



3 If you can. 
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FIGURE 3.8: Mass distribution resulting from a 3er fluctuation on 1 h^ 1 Mpc (solid 
lines). Black colour is used for e = 0.25, red for e = 0.5 and blue for pure 
radial orbits (e = 1). NFW (dashed lines) and Moore et al (dotted lines) 
profiles are shown for comparison. Left panel: Radial density profile. 
Right panel: Logarithmic slope qm of the mass profile. 



Moreover, since the radial coordinate of the particle changes only a fac- 
tor of three (|3.56j) . the approximation M(r) oc r Qju provides a reasonable 
description of the mass profile as long as the local value of the logarithmic 
derivative of the mass is used to set om(^)- 

In Figure ETHl the final profile obtained for a 3a peak in the primordial 
density field (smoothed on 1 h^ 1 Mpc scal es) is plotted for different v alues 
of the o r bit ecc entricity. As pointed out by White and Zaritskv ( 19921 ) and 
Hiotelisl ( 20fl2al ) . angular momentum plays a key role during secondary infall, 



preventing dark matter particles from reaching the centre of the halo and 
thus decreasing the amount of mass in the inner regions contributed by 
recently accreted shells. 

As the amount of angular momentum is increased (e is lowered), we find 
that the predicted density profile b ecomes increasingly shallow, in agreement 
with the results of lffiotelisl \2QQ2<h . Pure r adial orbits give rise to a profile 
somewhat similar to the form proposed by Moore et al. ( 1999bl ). although 
the exact shape is slightly different (see the right panel of Figure l3~%l where 
the logarithmic derivative of the cumulative mass is plotted as a function of 
radius). When the eccentricity is too low (as for example e = 0.25, repre- 
sented by solid black lines in the figure) , the halo develops a constant density 
core inconsistent with the results of numerical simulations. A value e = 0.5, 
typical of our clusters (see Section 13.2.2(1 , yields a mass distribution almost 
indistinguishable from the NFW formula up to 0.02 h^ 1 Mpc, where the 
logarithmic slope predicted by spherical collapse keeps increasing towards 
the centre while NFW tends to an asymptotic value au = 2. 

This is precisely the kind of behaviour we see in the innermost part of 
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our relaxed clusters (see Figure . We cannot be totally sure that this 
is not a spurious effect caused by our finite numerical resolution. There- 
fore, we suggest that a smoothing length of ~ 1 h^ 1 kpc would be needed 
to confidently distinguish between both profiles. We note, however, that 
this is only marginally beyond the resolution of the numerical experiments 
presented here. We plan to test this result with simulations from same ini- 
tial conditions but with another mass refinement level (i.e. effective 1024 3 
particles and e ~ 200j?c.) 



3.3.4 Comparison with simulations 



In principle, it is difficult to understand why the spherical collapse picture 
should be able to predict the density distribution of dark matter haloes, 
since the assembly of structures in the universe is not expected to proceed 
by the infall of spherical shells, but rather through a hierarchical series of 
merging events. Nothing of the symmetric and ordered process envisaged by 
this simple formalism is actually confirmed by the numerical simulations, in 
which the collapse of haloes seems to take place in an extremely unordered 
and random way. 



However, Zaroubi et al 



point out that the very complicated co- 
alescence process looks very regular in energy space, where there is a high 
correlation between the initial an final energies of individual particles. In 
particular, t he ranking in energy spac e is r oughly preserve d, an idea already 
proposed by Quinn an d Zurekj (Il988l'l and Hoffmaul ril98^ . 

Furthermore, iMoore et alJ (|1999bl ) simulated the collapse of a galaxy 
cluster truncating the power spectrum of primordial fluctuations at a scale 
~ 4 Mpc. The lack of small scale power causes the halo to form via a single 
monolithic collapse rather than through mergers and accretion of smaller 
haloes. Since the final density profile was very similar to that of the same 
halo with the standard power spectrum, they concluded that details of the 
merging history do not affect the final density profiles. 

Predictions of spherical infall model have been con fronted with n umer - 



ical simulation in several occasions. Early work by iQuinn et alJ (198 



found that in general the simple model corre ctly predicts t h e stru cture of 



Zarouhi et all £996) 



accom- 



bound objects, at least in a statistical sense, 
plished a series of 'simulations' in which forces were computed assuming 
spherical symmetry, obtaining the same results as a standard Tree N-body 
code (but, on the other hand, their analytic predictions did not agree s o 



well with the numerical resu lts). More recen tly, Lokas and Hoffman! (|2 



Del PopoIo et al 



200(1 1 and lHiotelisI (|2002al ) obtained density profiles com- 
parable to the NFW form using different implementations of the spherical 
collapse model. 

In order to test the validity of the formalism outlined above, we now 
try to reproduce the individual density profiles of our relaxed and minor 
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Cluster State u Rf 



A 
B 
C 
D 
E 
H 

Ji 

Ki 

L 

K 2 

■h 



Minor 
Minor 
Relaxed 
Minor 
Relaxed 
Relaxed 
Relaxed 
Relaxed 
Minor 
Minor 
Relaxed 



3.5 
3.4 
3.4 
3.2 
2.8 
2.5 
2.7 
2.7 
2.7 
2.5 
2.8 



1 
1 
1 

1 
1 
1 







3 
2 
3 
9 
4 
2 

9 
7 
9 
8 



Table 3.2: Spherical collapse parameters (peak height v and smoothing scale Rf). 

merging clusters. Assuming a fixed value of the orbit eccentricity, the model 
described in this section has two free parameters, which in this case have 
a clear physical meaning: namely, the smoothing scale Rf and height v (in 
units of sigma) of the primordial density peak from which the cluster is 
formed. 

In Table l3~2l we report the best-fit values found for our clusters. Most 
of them can be described as structures formed from ~ 3<r peaks of the 
primordial density field, smoothed on ~ 1 h~ l Mpc scales. It is interesting 
that there exists a certain degeneracy between the parameters Rf and v for a 
given mass. However, a higher peak on a smaller smoothed scale corresponds 
to an earlier formation time. The inner density profile of such peak would 
be much steeper than that of a similar object formed (at a later redshift) 
from a lower primordial peak on a larger smoothing scale. 

The average accuracy of the mass profiles predicted by our model based 
on spherical infall is shown in Figure 13.91 Individual profiles for each one 
of the clusters can be found in Appendix Although the overall quality 
of the fit is only slightly lower than that of the NFW fit (approximately 
10%, compared to 5 — 10% reported in Section l3.2.1|) . there exists a large 
scatter between different dark matter haloes. This additional uncertainty 
rises considerably the error in the mass profiles of individual clusters, which 
can be as large as ~ 25%. 

Nonetheless, the spherical collapse formalism is able to 'predict the mass 
distribution from a physically-motivated description, whereas the formulas 
proposed by NFW or Moore et al. are mere phenomenological fits (valid for 
a much narrower radial range). Therefore, it constitutes a promising tool 
to understand the origin of the so-called 'universal' density profile despite 
the fact that it does not include explicitly the merging effects on structure 



3.3. ON THE PHYSICAL ORIGIN OF DARK MATTER PROFILES 71 



0.1 1 
r [h _1 Mpc] 



10 



Figure 3.9: Accuracy of the spherical collapse prediction 



formation. 

We conclude that a possible explanation of why the density profiles pre- 
dicted by the spherical infall model might describe the real (or numerical) 
dark matter haloes could be that merging is taken into account implicitly 
through the smoothing scale Rf. Contrary to the common view (see e.g. 
BBKS), we argue that Rf has a very precise physical interpretation: dark 
matter haloes do not form in local maxima of the primordial density field. 
These are the seeds for the very first objects in the universe, which collide 
amongst themselves and give birth to heavier dark matter clumps. These 
clumps, in turn, correspond to peaks associated to a very small smooth- 
ing scale, below which all the information about the substructure of the 
primordial density field is lost. 

As time goes on, these protohaloes evolve through hierarchical clustering, 
and the smoothing length corresponding to the resulting objects increases 
(i.e. all the substructure below this scale is completely erased due to violent 
relaxation). Our point of view is that the smoothing length associated to 
a dark matter halo can be interpreted as a measure of the halo formation 
time or, equivalently, as the scale above which the system still remembers 
its initial conditions. 



CHAPTER 3. DARK MATTER 



Chapter 4 

Gas 



Too much pressure! Too much pressure! 
The pressure can't stop... 
Too much pressure! 

- The Selecter : Too Much Pressure (1980) - 



It the hierarchical galaxy formation scenario, cosmic structure grows 
via gravitational instabilities from small perturbations seeded in an early 
inflationary epoch. Since most of the mass is assumed to be in the form of 
collisionless cold dark matter, this yet unidentified component determines 
the dynamics of baryons on large scales, where hydrodynamic forces are 
unimportant compared to gravity. 

However, if we want our simulations to describe the observable properties 
of cosmological objects, we must obtain information about the physical state 
of the baryons, either by resorting to phenomenological relations between 
the dark and luminous components or by direct self-consistent numerical 
modelling. The most basic baryonic process to be included in a numerical 
simulation is gas hydrodynamics, i.e. the presence of a collisional term in 
the Boltzmann equation representing particle interactions. 

In principle, hot gas in clusters of galaxies should be easy to understand. 
Because of the relatively low ratio of baryons to dark matter, the potential 
of a cluster should be dark-matter dominated. The dynamical time within 
a cluster potential is shorter than a Hubble time, so most clusters should be 
relaxed. Also, the cooling time of vast majority of intracluster gas is longer 
than a Hubble time. It would appear that cluster structure ought to be 
scale-free, as long as the shape of a cluster's potential well does not depend 
systematically on its mass. 

If that were the case, then the global properties of clusters, such as halo 
mass, emissio n-weighted t emperature, and X-ray luminosity, would scale 
self-similarly (jKaiserl ll98fiy Indeed, numerical simulations that include adi- 
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abati c gasdynamics produce clusters of galaxies th at obey these scaling laws 
(e.g. Evrard et al. 199(4 Brvan and NormanI 1998h . 

Nonetheless, real clusters are not so simple. Even on cluster scales, ra- 
diative and feedback processes may not be completely negligible. This is 
expressed most clearly in the X-ray luminosity-temperature relation, prob- 
ably due to an increase in entropy of the innermost gas. Observations indi- 
cate a much steeper temperature dependence than that expect ed from pure 



David et al 



bremsstrahlung emission, especially for low mass systems (e.g 
ll 99,4 IPonman et ^1199^ . The amount of energy injection required de- 
pends on the gas density at the time when the heating occurred, either 
before, during or after the collapse. 

In this chapter, we address the reliability of the current numerical imple- 
mentations of gasdynamics, and then present our results for the sample of 
galaxy clusters described in Section T2.3.1I The ICM of dynamically relaxed 
clusters is found to be in hydr ostatic equilibrium and can be described by 
a polytropic equation of state ( Ascasibar et alJ 2002c ) . The radial profiles 
of gas density and temperature are also investigated, as well as the scal- 
ing relations between the global properties of our clusters. In most cases, 
the physical characteristics observed in the simulations show an encouraging 
agreement with relatively simple analytical predictions. 



4.1 Numerical gasdynamics 



Since the late 1980s a variety of techniques have been developed to simulate 
gasdynamics in a cosmological context. In part inspired by the success of 
the N-body scheme, the first gasdynamical techniques were based on a par- 
ticle representation of Lagrangian gas elements using the smoothed particl e 
hydrodynamics (SPH) technique ( Lucvlll977 : Gingold and Monag ianjjl97jtL 



S oon thereaf ter, fixed-mesh Eulerian methods were adapted (|Cen et al 



199 fli: ICenlll 992) and, mo re recently, Eulerian method s with sub-me s hing: 
dBrvan and NormanI 199.4 ) . deformable movin g mesh (lOnedinl 199.4 Pen 



1995, 1998) or adaptive mesh refinement (AMR. lBrvan et al .111994 : iKravtsov" et al 
2002 ) have been developed, as well as extensions of the SPH tech nique 
(|Shapiro et al~lll99fil : ISpringel and Hernauistll2002al : ISerna et al.ll2002h . 



These codes are actively being applied to a variety of cosmological prob- 
lems, ranging to the formation of individual galaxies and galaxy clusters to 
the evolution of Lya forest clouds and the large-scale galaxy distribution. 
Because the inherent complexity of gasdynamics in a cosmological context, 
such simulations are more difficult to validate than N-body experiments. 
Standard test cases with known analytical solutions (such as shock tubes) 
are far from the conditions prevailing in cosmological situations where the 
gas is coupled to dark matter, and this, in turn, evolves through a hierarchy 
of mergers. 
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The Santa Barbara cluster comparison project ( Frenk et al. 19991 ) at- 
tempted to assess the extent to which existing modelling techniques gave 
consistent and reproducible results in a realistic astrophysical application, 
simulating the formation of a galaxy cluster in a SCDM universe with 12 
completely independent codes, seven of them based on the SPH technique, 
while the other five employed grid methods to implement gasdynamics. 

The properties of the CDM component in all modes were encouragingly 
similar, most of the discrepancies being due to small differences in the timing 
of the final major merger, which happened to be uncomfortably close to 
z = 0. There was less agreement in the gas-related quantities, although 



in most cases they were relatively similar, usually within 10 — 20%. Only 
the X-ray luminosity showed a strong dependence on the resolution of the 
different codes, with a spread as high (IS (I factor of 10. 

The most remarkable finding was the hint of a systematic trend in the 
temperature profiles obtained for the inner regions (r < 100 kpc). Near the 
centre, SPH codes generated a flat or slightly declining (inwards) temper- 
ature profile, while grid codes produced temperature profiles that are still 
rising at the resolution limit. The entropy profile in SPH codes decreased 
continuously towards the centre, while grid codes develop an isentropic core 
at small radii. 



Several authors (e.g. Hernauistl 199.4 iNelson and Papaloizoul 199.4 19941 : 

3^ 



Serna et al . 1996; ISpringel and Hernauis t 2002 at ISerna et al 12 002 ) have pointed 
out that an important shortcoming of conventional SPH formulations is the 
poor conservation of entropy when a low number of particles is used. In 
this section we intend to test the accuracy of both Eulerian and Lagrangian 
numerical approaches for integrating the equations of gasdy namics. To this 



end, we compare results of the adaptive mesh code ART (jKravtsov et al 



3, (see Section HOI for a description) with two different implement, 
tions of the Lagran gian code Gadget; on e of them is the publicly-available 
version of the code (jSpringel et al. 2001bl ) , which implements standard SPH 
gasdynam ics, and the other is th e explic it entropy-conserving scheme de- 
scribed in Springel and Hernquist ( 2002al . see also Section 12. Q[) . 



4.1.1 Santa Barbara cluster 

As a first step to test our codes, we have chosen the Santa Barbara clus- 
ter initial conditions, corresponding to a 3a peak in the primordial density 
field smoothed on a 10 Mpc scale. Details concerning the numerical exper- 
iments accomplished at different mass and spatial resolutions can be found 
in Section T2.3.2I 

We have plotted in Figure 1-01 the results of all these experiments, using 
the criterion of 100 gas particles, 200 dark matter particles and r > 3e to de- 
fine the minimum radius i n each case. In partic ular, we compare the results 
of the Eulerian code ART ( Kravtsov et al.ll2002l . triangles) with the standard 
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FIGURE 4.1: Simulations of the Santa Barbara Cluster. Top left: Dark matter den- 
sity. Bottom left: Gas density. Top right: Gas temperature. Bottom 
right: Gas entrop y. Solid squa res denote the average of the Santa Bar- 
bara experiment l|Frenk et a 1 .1 Il999h . Lines represent Gadget results 
for simulations with different mass resoluti ons and smoothing (see text), 
triangles are used for ART and circles for (|Brvan et alJ ll995). 



and entropy-conserving implementations of the SPH algorithm available in 
Gadget (re d and blue lines, re spectively). The average of the simulations 
presented in lFrenk et al.l (|1999l ) is sho wn by the sol i d squ ares. The results 
obtained by the adaptive-mesh code of iBrvan et al. (1995) are plotted sep- 
arately (circles). 



Resolution 



The profiles for the different numerical experiments of the Santa Barbara 
cluster reported in Table l2~31 can be seen in Figure l4~Tl Blue lines are used 
to represent the simulations with 256 3 particles and e = 2 h~ x kpc (dotted), 
128 3 - 5 h- 1 kpc (dotted-dashed), and 128 3 - 20 h' 1 kpc (solid). Red 
lines (standard SPH) correspond to 128 3 - 20 hr 1 kpc (solid), 128 3 - 5 hr x 
kpc (dotted-dashed), and 64 3 — 20 hr 1 kpc (dotted). The black line on 
the bottom right panel indicates the entropy profile obtained in the low 
mass resolution simulation with 64 3 particles run with the standard SPH 
implementation of Gadget, ignoring the resolution limit. 

In most cases, all the res ults are consi s tent when plotted according to 
the restrictions prescribed by Klvpin et alJ (|200ll ). Nonetheless, the two ex- 
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periments with the highest nominal spatial resolution (5 hr 1 kpc and 256 3 
particles with standard and entropy-SPH) develop a steep decrease in the 
central temperature. Since we have checked that the centre of mass of both 
gas and dark matter distributions are approximately coincident (see e.g. the 
high gas densities obtained for the standard SPH implementation), this ef- 
fect is probably caused by an artificially low smoothing length, which leads 
to very compact groups of cold particles that decouple from the surrounding 
medium. Although Eulerian simulations are not affected by this problem, 
the unphysical temperature drop described here is often seen in cl usters sim- 
ulated with SPH-based code s at low and medium resolutio ns (e.g. lEke et al. 
1998; lYoshikawa et all 1200(1 Mathiesen and EvrardlBoOll ). and in some of 
the codes reported in lFrenk et alJ (|l99*9 ) . 

If the smoothing scale is not set excessively low, but a small number of 
particles is used, the resulting cluster displays an extended core of constant 
density and temperature up to the minimum radius we have used for our 
plot. This softening-dominated core yields a flat entropy profile at the inner 
regions (black line in Figure l4~Tj) . pointing to the misleading conclusion that 
this feature of clusters and groups of galaxies can be reproduced with a fairly 
low resolution. The s ame effect has been s hown in cluster simulations with 
variable resolution by Borgani et al. ( 2002h . 



Code comparison 



A different issue concerning the entropy profile of galaxy clusters is whether 
the computations based on different integration schemes give consistent re- 
sults, for th e same numerical resolution. This topic has been previously 



discuss ed by Kan g et al 



( 20023) or lSerna et al 



( 1994) , iFrenk et al.l Jl999) , S pringel and Hernquisli 
2002), among others. 



The maximum number of particles allowed by the Santa Barbara initial 
conditions was 256 3 , which seems insufficient to discern the presence of a 
physical flattening of the entropy from the effects of nume rical smoothing. 
However, we confirm the trends shown in lFrenk et alJ (|1999i ): Eulerian codes 
predict a systematically rising temperature profile near the centre, while 
SPH algorithms are more consistent with isothermality. 

This conclusion is strengthened by the fact that the two completely 
independent grid- based methods used in this comparison (i.e. ART and 
Brvan et al. 199fj ) give very similar results, while both imple mentations of 
Gadg et show profiles consistent with the average value given in lFrenk et al. 
(|l999h . which was dominated by SPH experiments. 

Nevertheless, it seems that the highest resolution simulations using the 
entropy-conserving version of Gadget develop a certain tendency towards 
higher central temperatures, displaying an increasing profile up to ~ 60 h~ x 
kpc. At this point, there is a sharp break and the temperature drops again in 
the inner part of the cluster. Quite interestingly, the entropy profile at that 
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Figure 4.2: Simulations of Cluster A 



radiu s shows remarkable signs of flattening (a fact also noted bv lSerna et al 



2002, who try to overcome the entropy problem by including the so-called 
V/i terms arising from the spatial dependence of the SPH smoothing length). 
In any case, more resolution is needed in order to study the structure of the 
ICM gas at even lower scales and make a reliable assessment on whether the 
entropy profile produced by these improved formulations of SPH is consistent 
with the results of Eulerian codes. 



4.1.2 Cluster A 

Following the spirit of the Santa Barbara Cluster Comparison Project, we 
have run several numerical experiments from the initial conditions that pro- 
duced the most massive galaxy cluster of our catalogue, as described in 
Section In this case, the number of particles within the virial radius is 
approximately a factor of 8 larger than in the Santa Barbara cluster. Hence 
we have been able to resolve scales much deeper into the cluster's potential. 
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Figure 14.21 shows the results obtained with ART (|Naeai and Kravtsov 



2002) and both versions of Gadget for the radial profiles of Cluster A. 
Here, the difference between energy and entropy conservation is mostly ev- 
ident at low radii (below ~ 20 h^ 1 kpc). As noted in previous studies (e.g. 
berna < ^HD[2002), the density profile is steeper in the standard formulation 
of SPH, whereas the temperature is systematically higher when entropy con- 
servation in SPH is enforced. 

Consequently, the entropy profile is dramatically dependent on the cho- 
sen implementation. In top right panel of Figure 14. 2[ it can be clearly 
seen that the entropy profile o f the explicit entropy-conserving scheme by 
Springel and HernquistJ ( 2002al ) unambiguously flattens in the central re- 



gions of the simulated cluster, while the standard SPH decreases uninter- 
ruptedly down to the resolution limit. 

Since this effect is already manifest at clustercentric distances that are 
confidently resolved, we claim that there is a fundamental difference in the 
predictions based on Eulerian and Lagrangian codes due to the unphysical 
entropy loss in the latter. The fact that the entropy-conserving version of 
Gadget produces an entropy floor comparable to that observed in ART 
gives reassuring support to the idea that this numerical artifact is corrected 
in this formulation of the SPH algorithm. 



4.2 Structure of the ICM gas 

It has been twenty years since the reali zation that the extended X-ray emis- 
sion from clusters ( Kellogg et alj 1972h is thermal a nd arises from optically 



thin plasma filling the clusters ^Mitchell et ah 197Gh . If this plasma is in hy 



drostatic equilibrium within the gravitational potential created by the dark 
matter, the spatial distribution of the ICM density and temperature con- 
stitutes an invaluable source of information about the process of structure 
formation on large scales. 

Most importantly, since clusters of galaxies are the largest gravitationally 
bound structures in the universe, their physical properties have profound 
cosmological implications, and have often been proposed as a probe to test 
the values of the cosmological densit y parameters as w ell as the validity of 
the hierarchical CDM scenario (e.g. White et al. 199.il ). If well calibrated, 



the slope and evolution of cluster scaling relations, such as size, gas mass or 
luminosity versus temperature, can be used to constrain cosmological and 
structure formation models. Temperature profiles are also a fundamental 
ingredient in the determination of the total mass, as well as the gas entropy 
distribution, which is a powerful tool to explore non- gravitational processes 
that could alter the specific thermal energy in the ICM. 

It is therefore important to investigate in detail the structure of the 
intra-cluster medium. Throughout this section, we will assume spherical 
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symmetry and focus on the radial profiles of several quantities of physical 
interest: first, we will address the validity of the usual assumptions of hy- 
drostatic equilibrium, isothermality, and polytropic equation of state. Then, 
the average density and temperature profiles of our galaxy cluster catalogue 
(whose main properties are described in Section r2J-i.il) will be compared with 
predictions based on the existence of a universal density profile for the dark 
matter component. 



4.2.1 Physical state 

One of the hypothesis that are often used in the study of galaxy clusters is 
that gas is in hydro static equilib rium with the mass distribution dominated 
by the dark matter. Alien! ( 1998h suggested that this assumption is only valid 



for cooling-flow clusters, which has been corroborated by recent observations 
with Chandra and XMM-Newton. 

Most analytical studies on the radial profiles of the ICM gas rely on 
hydrostatic equilibrium, in many occasions supplemented by a polytropic 
or isothermal equation of state for the gas. As we briefly discussed in Sec- 
tion EH1 these prescriptions are used, for example, to estimate the grav- 
itating mass from the observed X-ray emission. It is therefore of crucial 
importance that the accuracy of these assumptions is investigated by means 
of direct numerical simulation. 

The main difficulty in studying the properties of gas in clusters formed in 
nume rical experi ments is the of lack o f resolution (e.g. lThacker and Couchmanl 
200l]). Recently, lLoken et al.l twA ) have made a detailed analysis of the 



gas temperature profile for a sample of 20 clusters simulated with a resolu- 
tion (< \5ti~ 1 kpc) which is comparable to the present work, using an AMR 
gasdynamical code. 

A completely different problem is the implementation of all the relevant 
physics involved in cluster formation. Radiative cooling and stellar feedback 
play a major role on galactic scales, as well as in the dense central parts of 
galaxy clusters. Outside the core, though, cooling times are longer than a 
Hubble time, and the effect of cooling on the ICM structure is not expected 
to be significant, at least on a first approximation. 

On the other hand, there are still many open questions regarding the 
assembly of galaxy clusters even when only adiabatic physics is taken into 
account (see e.g. the previous section). Before including any radiative pro- 
cess, it is necessary to make sure that the physics associated to gravothermal 
collapse and shock- wave heating are accurately modelled. This also applies 
to the many theoretical predictions that can be derived with this 'simple' 
physics, which are often used as a benchmark to be compared with the 
observed properties of real clusters of galaxies. 
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Figure 4.3: Test of the hydrostatic equilibrium assumption 



Hydrostatic equilibrium 

Assuming that the gravitational potential is counter-balanced by thermal 
support only (i.e. bulk - infall and/or rotation - and turbulent motions of 
the gas are negligible with respect to its internal energy), the hydrostatic 
equilibrium equation for a spherically symmetric system reduces to the sim- 
ple form 

J_dP _ _CMM 
Pg dr r z 

where p g (r) is the gas density, P(r) = p3 ^ kT ^ is the pressure and M(r) is 
the total (gas + dark matter) mass enclosed within radius r. Since the dark 
component dominates the mass of the cluster at all radii, the cumulative 
mass can be approximated by any of the 'universal' profiles described in 
Section 13.1.11 

We have checked the reliability of the usual assumption of hydrostatic 
equilibrium in Figure 14, 3( where the ratio of the pressure gradient to the 
gravitational term (pGMr~ 2 ) is plotted as a function of clustercentric radius. 
If the gas was in thermally-supported hydrostatic equilibrium, this quotient 
should be equal to unity by virtue of equation 1)4. ip . 

For clusters that have been classified as Relaxed systems on dynamical 
grounds (see Section I2.3.1|) . hydrostatic equilibrium can be applied with 
a reasonable degree of accuracy (~ 10%). Deviations from the behaviour 
predicted by (|4.1|) are due to the oversimplifying assumption that angular 
momentum and turbulent motions do not contribute significantly to support 
the gas against gravity. Besides, near the virial radius (r > 0.8r2oo) 5 radial 
infall cannot be completely neglected. 

As can be seen in the middle and right panels of Figure 14.31 the assump- 
tion of hydrostatic equilibrium holds only marginally for minor mergers, and 
not at all for clusters that are undergoing a major merging event. At this 
point, it is important to recall that non-relaxed systems constitute 60% of 
our sample. Even if the error is not extremely large (particularly, in the case 
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Figure 4.4: Polytropic index as a function of r 



of minor mergers), caution must be kept in mind whenever equation (|4.1|) is 
applied to a real cluster ensemble. 



Polytropic equation of state 

In many occasions, hydrostatic equilibrium is used in combination with a 
simple equation of state. Usually, the gas is assumed to follow a polytropic 
law 

p = PjM_ = 

firrip y 

where 7 is the polytropic index and the normalisation Pq is an arbitrary 
constant. The value 7 = 1 corresponds to the special case of an isothermal 
gas, and is often seen in the literature, both in analytical and observational 
studies. 

If the ICM gas could be described by the polytropic form ()4.2j) . the radial 
temperature and the density profiles are related by 



To \po 



(4.3) 



where po and To are the central values of both magnitudes. 

In real clusters, there is no strong physical reason to expect that the 
gas should be described by a polytropic equation of state. However, the 
appealing simplicity of equation (|4.2|) makes it extremely useful in analytical 
work, and hence it would be interesting to test whether it can be used to 
interpret the results on a first-order approximation. 

The most straightforward way to check whether the ICM follows a poly- 
tropic relation is to compute the local value of the polytropic index 



(4.4) 
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Figure 03 displays the radial dependence of 7 for our sample of simulated 
galaxy clusters. Albeit some scatter, the radial profiles of the diffuse gas 
component in relaxed systems and minor mergers seem to be consistent 
with a polytropic equation of state with 7 ~ 1.18 (shown by a horizontal 
dotted line in the figure). 

Another important issue is that an isothermal profile (7 = 1) offers a 
very poor approximation to our results. In the section below, where the 
temperature structure of our clusters will be described in more detail, we 
will see that the constant temperature approximation can be ruled out at a 
significant level of several sigmas. 



4.2.2 Radial profiles 



In the previous chapter, we thoroughly discussed the possibility that the 
radial mass distribution of dark matter haloes followed a 'universal' form, 
finding that relaxed haloes (and even minor mergers, albeit with a lower 
degree of accuracy) could be described by the simple analytical fits proposed 
in the literature. 

Once we know the dark matter distribution, the hydrostatic equilibrium 
equation (|4,1|) relates the gas density and temperature profiles to the cumu- 
lative mass, which can be approximated by the dark matter mass. 

Several further assumptions must be made in order to compute the gas 
and temperature profiles. By f ar, the most common prescrip t ion is the so- 
called /3-model, introduced by Cavaliere and Fusco-Femianol ()l97fil ) to de- 
scribe the X-ray surface brightness profiles of galaxy clusters. Assuming 
that gas is isothermal and in hydrostatic equili brium with a dark matter 
potential that follows an analytical King profile (|Kinglll962l )). 



p{r) 



Po 



r 

{1 + x 2 ) 3 / 2 ' X = 7 C 



M(r) = A-Kporl 

(j)(r) = -AkGpqt^ 
one can derive the gas density profile 

Pg( r ) 



ln(x + yl + x 2 ) 



x 



ln[x + (l+x 2 ) 1 /2] 



[l + (r/r c ) 2 ] 3 ^ 2 ' 



(4.5) 
(4.6) 
(4.7) 

(4.8) 



where po is the central density, r c is a core radius and (3 ~ kT / {pm p a 2 DM ) 
relates the thermal energy of the gas to the velocity dispersion of the dark 
matter component. This form for the density profile p g has been widely used 



in the literature (see e.g. iRosati et al.ll2002l .for a review). More recently, a 



generalisation of the /3-model to account for temperature gradients in the 
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surfac e brightness profile of a polytropic gas has been derived by Ettoril 



As we have shown in the previous chapter, the King profile is not the 
model that best fits the cluster gravitational potential of dark matter halos 
in numerical simulations. 



Suto et al.l (J1998) applied the hydrostatic equilibrium equation to a va- 
riety of dark matter density distributions. Assuming a polytropic equation 
of state for the gas and a NFW profile for the dark component, they find 

m^n^^pTj ( , 9) 

and a similar (though more complicated) expression for a Moore et al. for- 
mula. The B parameter gives the normalisation of the temperature profile, 
To, as a function of cluster dark mass, parametrised in terms of the charac- 
teristic density p s and radius r s of the NFW profile: 

7 k B T 

Needless to say, the corresponding gas density profiles are given by (|4,3|) . 
B fixes not only the central value To, but also the asymptotic behaviour of 
both gas and temperature profiles (T(oo) — ► 1 — B). If we want the baryon 
fraction to remain approximately constant at large radii, we must choose 
B = 1 and 7 ~ 1.2 in order that p g ~ r~ 3 . 

A similar result for B and 7 has been found by Komatsu and Seliakl 



(2001) from the same assumptions (i.e. NFW profiles, hydrostatic equi- 
librium, polytropic equation of state and constant baryon fraction at large 
radii). The normalisation of their temperature profiles is defined as 

3^To^ = 3 7_ r _lc^) 
/v ; GM 200 pm p 7 B y J 

where g(c) = [ln(l + c) — j-^] . However, these authors did not require 
that B = 1 to have a constant baryon fraction. Instead, they propose a 
relation between the normalisation r?(0), the polytropic index 7 and the 
concentration parameter c, which can be fitted by 

^(O) = 0.00676c 2 + 0.206c + 2.48 (4.12) 

and 

7 =1.15 + 0.01c (4.13) 

wi th c = cnfw — 6.5. T h is fit w orks well for c > 3, as can be seen in Figure 3 
of Komatsu and Seliakl ( 2001 ) , but cannot reproduce the normalisation for 



less concentrated halos. 
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Figure 4.5: Left panel: Temperature normalisation ?j(0). Right panel: Concentration 
dependence of the polytropic index. In both cases, solid lines are used for 
1[4.14H and dashed lines for 7 = 1.18. Dotted lines represent the results 
of Ko matsu and Selial3 | 



In Figure l4"31 we compare Q4.12JI with our prescription B = 1, equivalent 
to 77(0) = 3^^-cg(c). Assuming a constant polytropic index, 7 ~ 1.18, as our 
simulations seem to indicate (see Figure UU, we obtain a slightly different 
normalisation. In order to obtain a similar prediction for r/(0), we should 
use the following dependence of 7 on the concentration parameter: 

7 = 1.172 + 0.006c (4.14) 

As can be seen in the figure, the match with the analytical fit (|4.12j) is 
perfect fo r c > 3. For lower values of c , our prescription follows the numerical 
results of lKomatsu and S elia (see Figure 3 of their paper) . 



In any case, both models are compatible with results from our simulated 
clusters. The most important conclusion drawn from these analyses is that 
physical properties of the gas are intimately related with the dark matter 
distribution. As we will show below, the radial temperature and density 
profiles of the ICM can be unambiguously determined from the parameters 
of the dark matter halo. Therefore, it is expected that they can also be 
characterised by 'universal' formulae. 

Temperature 

Mass determinations from X-ray emission in clusters usually assume that 
relaxed clusters (i.e. morphologically symmetric) are isothermal. However , 
this assumption has been questioned recently bv iMarkevitch et al. (1998), 



who found evidence of a decreasing temperature profile for nearby clusters 
observed with ASCA. A large subset of their sample showed a remarkably 
self-similar temperature profile when properly normalised and rescaled to the 
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Figure 4.6: Average temperature profile (dots) compared to the theoretical predic- 
tions based on NFW (red line) and Moore et al. (blue line) profiles for the 
dark matter, assuming hydrostatic equilibrium and a polytropic equation 
of state. Accuracy of the NFW-based fit is shown in the top panel. 



virial radius. On the contrary, White! ( 20001 ) and Irwin and Bregman (200 



using data from BeppoSAX and ASCA satellites, do not find any decrease of 
the temperature in a large collection of clus ters. More recently, the analysis 
of BeppoSAX observations accomplished bv lDe Grandi and Molendil ( 2002) 



concluded that temperature profiles of galaxy clusters can be described by 
an isothermal core followed by a rapid decrease. 

In this context, we would like to investigate what can we learn from the 
temperature profiles obtained in our numerical simulations. Particularly, 
we will focus on the shape of the radial temperature profiles, compared 
with both observational data and the simple analytical predictions based on 
hydrostatic equilibrium and a polytropic equation of state. 

The mass-weighted temperature of our sample of galaxy clusters (ex- 
cluding major mergers) is shown in Figure H~B1 as a function of radius. We 
see that, with this normalisation, the scatter is low enough (~ 15%) to hint 
that the temperature profile can be described by a 'universal' form. 

The scale factors in both axis of Figure 14.61 have been chosen to com- 
pare with the predictions based on a NFW profile, hydrostatic equilibrium, 
polytropic equation of state and approximately constant baryon fraction at 
large radii (i.e. B = 1, T(oo) — > 0): 

T(r) ln(l + — 1 v - 1 

-P-= — rJ ! k B T = l ^Gixm vPs rl (4.15) 

T - 7 
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Figure 4.7: Averaged projected temperature profile (b lack squares with error bars) 
compared to the observational results of lDe Grandi and Molendi 
(solid circles for cooling-flow clusters, empty cicles for non-cooling flow) 
and iMarkevitch et al.l lll99 ? l) (dashed blue line). Red lines show the 
analytical fit proposed bv lLoken et al.l (|2002h (see text). 



which is shown as a red line in Figure l-PA In the top panel, we plot deviations 
from this 'universal' profile. The temperature profile derived from Moore et 
al. density profile is also shown for comparison (blue line). As can be easily 
appreciated in the figure, the difference is very small even near the centre. 

We would like to remark that equation Q4.15JI is not a fit to the radial 
temperature profiles of our clusters, but a prediction based on the existence 
of a universal dark matter distribution (in this case, given by NFW ex- 
pression). The central temperatures To have been computed assuming a 
polytropic index 7 = 1.18 for every cluster, so there are no free parameters 
related to the gas distribution. When only relaxed clusters are considered, 
the scatter around the predicted profile (determined by the values of p s and 
r s ) is considerably reduced. 

In spite of its physical significance, the radial temperature profile has 
the disadvantage of not being a direct measurable quantity. In order to 
make a more reliable comparison with observations, the projected emission- 
weighted temperature profile s must be computed from ou r numerical data. 
More elaborate methods (e.g. Mathiesen and Evrardl 200l] ) could be used to 
mimic the instrument-dependent observational procedure, simulating annu- 
lar spectra and estimating the temperature profiles in the same way that is 
done in real clusters. 
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Recent resu lts from high reso lution AMR gasdynamical simulations of 
galaxy clusters ( Loken et al.ll2002l ) seem to indicate that the projected X-ray 
temperature profiles have indeed an universal form which decreases signifi- 
cantly with radius. These authors propose, as a best fit to their numerical 
simulations, the simple analytical formula 



T(r) 



(1 + rKf 



(4.16) 



where To is the central temperature and a x is a core radius. In order to 
confront this exciting result with our numerical experiments, we compare, 
Figure l4~Tl the averaged temperature profile (once again, excluding clusters 
with major mergers) projected along the line of sight normalised to the 
averaged X-ray temperature at the virial radius 1 . 

Solid red line shows the best fit we obtain for our data using expression 
(14.161). wher e as the dotted red line indicates the radial dependence found by 
Lok en et al.l (|2002l) . Although the results seem to be consistent within the 
error bars, and the normalisation is similar in both cases (To — 1-33 (Tx)), 
the best-fit values of the core radius and the asymptotic expo nent are ex- 
treme ly sensitive to the particular details of the profile. While lLoken et al.l 
(2002) quote a x = r2oo/l-5 and 5 = 1.6 in their study, our clusters are better 
fit by a x = r2oo/4.5 and 5 = 0.7. 

As can be seen in Figur e 11*771 our results ar e almost identical to the obser- 



vational estimates of lMarkevitch et al.l ([19981 ) . who clai m that they a re we ll 



represented by a polytropic /3-model. But, as noted bv lLoken et al. (2002), 
this is not entirely self-consistent, since gas in the /?-model is assumed to 
be isothermal and the data show a pronounce temperature gradient. Nev- 
ertheless, because gas in our cluster do follow a polytropic relation between 
density and temperature, and we expect the density of gas to be close to the 
observed one, it is not surprising that the projected profiles agree so well 
with observations. 

Howe ver, the agreement is not so goo d with the temperature profile ob- 
tained bv foe Grandi and Molendl (|2002l 'l for a selection of clusters observed 
with BeppoSax. Our clusters are marginally consistent with their data, but 
we do not find any evidence of the temperature flattening observed by these 
authors at the inner regions. This could reflect the need of additional non- 
adiabatic physics in the simulations, but it could also be a consequence of 
the observational difficulties involved in the spectral determinatio n of the 
gas temperature (see the discussion in Mathiesen and Evrardl 2001 ). 



1 For a detailed description on the computation of these quantities, the reader is referred 
to Chapter "3 
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Density 

Assuming that the ICM gas can be described by a polytrope, the gas density 
profile in hydrostatic equilibrium with a NFW dark matter distribution is 
given by the formula 



Pg( r ) 
Po 



'ln(l + D 



7-1 



(4.17) 



where we have taken B = 1, as we did in (|4.15|) for the temperature. With 
this prescription, the gas density vanishes as r — > oo and the logarithmic 
slope is 

1 



a(r) 



7-1 



-1 + 



(4.18) 



-1 + 



Given the asymptotic behaviour at large radii, a ~ 

it is impossible that the gas density does exactly follow the dark matter 
pr ofile in order to maintain a constant baryon fraction (a fact already noted 
bv lKomatsu and Selia However, we should recall that NFW is not 

intended (nor able) to fit the dark matter density much beyond the virial 
radius, and the assumption of hydrostatic equilibrium is no longer valid in 
the outer parts of the cluster (see Figure FO]) . Therefore, the only, and 
necessary, requirement for (|4.17|) to be a good approximation to the ICM 
density is that it had a similar shape to the NFW profile for a broad radial 
range arou nd r^nn- 



This led lKomatsu and Seliakl ( 200ll ) to the constraint (|4.13j) on the poly- 
tropic index of the gas, but no attempt was made to obtain a density normal- 
isation, arguing that the gas mass fraction does not appear in the hydrostatic 
equilibrium equation, and that t here is observational evide nce suggesting its 
dependence on cluster mass (e.g. Arnaud and Evrardl 1999h . This makes the 
gas fraction determination uncertain and hence indicates that it should be 
treated as an additional free parameter. 

Although the global baryon fraction is not a constant among the different 
clusters, we propose that a first-guess approximation to the central density 
po may come from the local ratio of the gas and dark matter densities. 

This quantity is plotted in Figure 14.81 as a function of clustercentric ra- 
dius, normalised to the cosmic baryon fraction f2b/^dm- Not surprisingly, 
the radial profiles of both dark and gaseous components are very similar 
throughout most of the cluster 2 , differing only near the origin, where pres- 
sure is responsible for the formation of a shallow gas core. A flat baryon 



fraction profile has also been observed by Allen et al. (2002), who applied 



this result to the estimation of the mean matter density of the universe. 

In fact, this was one of the basic assumptions in the work of iKomatsu and Seliakl 
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Figure 4.8: Radial dependence of the local baryon fraction 



The baryon fraction expected from our prescription (|4,17|) . plotted as a 
red line in Figure l4~8l is simply 



Pg\ 



i 

7-1 



Pdm(t) 



-(1 + -) 



(4.19) 



Unfortunately, this expression does not tend to a definite asymptotic 
value, but it varies slowly enough at large radii for the range of polytropic 
indexes in which we are interested (as 7 = 1.18 in the figure). Actually, the 
shape of the gas density profile would be better modelled by slightly higher 
values (7 ~ 1-2), resulting in a more constant baryon fraction in the outer 
parts of the clusters, but we will continue to use same value of 7 as in the 
temperature profiles in order to maintain self-consistency. 

The averaged radial density profile of the ICM gas is shown in Figure l-Ol 
rescaled by r s on one axis and the value of po inferred from (|4,19|) in the 
other. The normalisation po has been obtained by imposing the condition 
that Pg/pdm = Qb/^DM at the radius where (|4.19|) has a maximum. For 
7 = 1.18, that is 

Po - 1.514^/), (4.20) 



dm 



The quality of this approximation, given in the top panel for the density 
profile derived from NFW (red line) is significantly poorer (~ 20 — 30% 
accuracy) than the phenomenological fits for the dark matter. In this plot, 
we also show the gas density expected from a Moore et al. profile (blue line). 
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Figure 4.9: Average gas density profile 



As in the temperature, the difference with NFW is not very large (usually, 
of the same order as the scatter around the NFW-based prediction). Part 
of this scatter can be attributed to departures from hydrostatic equilibrium 
in merging clusters, as well as deviations from purely thermal support that 
take place even in the outer parts of relaxed clusters. 

Apart from self-consistency, another advantage of parametrisation Q4.17|) 
over the conventional /5-model is due to the instability of the parameter 
with respect to changes in the outer radius used to infer its best-fit value. 
Observational data restricted to radii much smaller than r2oo tend to find 
values close to (3 ~ 2/3, wh ile wider field cluste r images are usually better 
described by f3 



1 (see e.g. Navarro et al. 19951 ). Since the profile given in 



expression (|4.lYj) depends only on the dark matter parameters, and these 
are much more stable than those of the /3-model, we consider this formula 
to be more reliable in order to extrapolate the gas density up to the virial 
radius. Moreover, it allows a direct estimate of the total mass if the dark 
matter halo can be described by a NFW profile, although we have already 
shown that the accuracy is somewhat limited. 



4.3 Scaling relations 

Before the structure of clusters of galaxies could be resolved by the available 
observations, they were already realised to exhibit a number of simple cor- 
relations between their global properties. These scaling relations could lend 
insight into the physical nature of clusters when compared with the analyt- 
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ical predictions expected to hold under the assumpti on of self-sim ilarity in 
the formation and evolution of large scale structure 

However, deviatio ns from self-similarity are expected under the effects of, 
for example, merging (.Ti ng and Su to 2000) and any additional physics acting 
on t he intracluster gas oyer the simplistic infall in the dark matter potential 
(e.g. Evrard and Henrvlll991 : brvan and Normanl 19981 : Bialek et al. 200ll . 
and references therein). The latter case is particularly relevant to cool sys- 
tems, where the extra energy es timated from observed properties is compara- 
ble to their thermal energy (e.g. Ponman et al.ll99fil . 19991 : iTozzi and Normanl 
200lh . 

In the past 10 years, the advent of space-based X-ray observations of 
galaxy clusters have revealed that there exist tight correlations between 
the total gravitating mass, th e X-ray luminosity and the temperature of 



the intracluster medium (e.g. Rosati et al. 2002). A close comparison of 



theoretical results and observational data can therefore provide valuable 
constraints on the prevailing cosmological models and even on the nature 
of the dark matter that dominates the mass distribution and the dynamical 
evolution of clusters. 



4.3.1 Mass-temperature 

The number density of galaxy clusters as a function of their mass can provide 
a powerful probe of models of large-scale structure. Since the mass of a 
cluster is not a directly observable quantity, the abundance of rich clusters 
of galaxies has been historically measured in terms of some other parameter 
which is used as a approximation for mass. Several options exist, but much 
attention has been focused recently on the X-ray temperature. 

Cosmological simulations suggest that this magnitude is strongly corre- 

lated with mass, showing little scatter (e.g. Evrard et alll99"fil : lBrvan and Normanl 
1998). How well simulations agree with observational results is still far from 
clear, and several issues need to be resolved. Numerical resolution and 
difficulties to include all the relevant physics are the main sources of un- 
certainty invariably coupled to cosmological simulations, while instrumental 
effects and the lack of a reliable method to estimate the mass are the most 
worrying aspects related to the observational determination of the M — T 
relation. 



Analytical prediction 

Within the framework of pure gravitational infall (e.g. lLilielll992l ). the clus- 
ter mass and the gas temperature should scale simply as M oc T a , where the 
exponent a ~ 3/2 is almost unrelated to the particular settings of cosmo- 
logical paramete rs. This self-similar i ty has been later justified by numerical 
simulations (e.g. lEvrard et a~ 199fil : Brvan and Normanl 1998) and is often 
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used as a reference against which observational results are compared. 

There are several ways to obtain this relation. The simplest of them 
relates the thermal energy of the gas to the potential energy due to self- 
gravitation of the cluster. By virtue of the virial theorem: 



kT A 
fj,m p 



I GM A 
2 r A 



(4.21) 



where T and M refer to average temperature and cumulative mass. The 
subscript A indicates the overdensity threshold at which these quantities 



are measured. Substituting M A = ^-Ap c r A and p c = |^ above we get 



M 4 fkT A \ 3/2 

A ~ GHqA 1 / 2 \pm p ) 



1 keVy 



(4.22) 



Therefore, the virial theorem predicts not only the relation between M and 



T 1,5 , but also the normalisation M { 



vir 




A similar scaling law can be derived from our assumption of a poly- 
tropic equation of state for the gas, in hydrostatic equilibrium with a NFW 
potential. Using the density profile (|4.17|) . the averaged temperature at 
overdensity A is just 



T A = T T m (j,r A /r s ) ; T m ('j,r A /r s ) = 
and the emission-weighted temperature 
Tx = T T x (j,r A /r s ) ; T x (>y,r A /r s ) 
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(4.24) 

The assumption B = 1 led to the central temperature normalisation 
given in equation (|4.15j) . Using the definition of M A , we recover the scaling 
relation M A = MqA' 1 / 2 ^ 2 , with M given by 



M n 



Mq 11 



1 



ln(l + r A /r s ) 



TA/r s 
I+ta/t 



r A /r s 



3/2 



(4.25) 



where the term in brackets depends on the particular values of the polytropic 
index and 'concentration' r A /r s of each dark matter halo. As we discussed 
in the previous section, the polytropic index can be taken to be a constant 
value 7 ~ 1.18 for the range of masses covered by our simulated clusters. 
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Figure 4.10: Analytical M — T normalisation. Top panel: Mass-weighted tempera- 
ture. Bottom panel: Emission weighted temperature Tx- 



The validity of this approximation outside this mass range must be tested 
by independent numerical simulations. 

The dependence of Mq on r^jr s is shown in Figure E.l()[ which can be 
interpreted as a 'universal profile' of the ratio between thermal and grav- 
itational energy. The normalisation Mq = MaA x / 2 T a 3//2 is plotted as a 
function of radius, in units of the characteristic scale r s . For a value of the 
NFW concentration parameter c = 8, the radius at which the overdensity is 
A = 2500 approximately coincides with r s /3, riooo ~ r s/2 and rsoo ~ 2r s /3. 
By definition, r2oo = cr s . 

As can be seen in Figure 14. the normalisation of the M — T rela- 
tion predicted by equation (|4,25|) depends on whether T or Tx is used, the 
difference (which can be up to a factor of 2 near the virial radius) is due 
to the factor t m x . Our formula gives a fair approximation to the mass- 
temperature relation in relaxed clusters, but it severely overestimates the 
mean temperature in merging systems (by > 50% in the outer parts). At 
high over densities, though, expression (|4.25|) is consistent with the results of 
our simulations within the scatter from halo to halo. The effect of dark halo 
concentration (which varies between 7 and 10 for our relaxed clusters, 5 and 
7 for minor mergers, and 3 — 5 for the major merging subset, see Table l3~T|) 
is always of the order of 10% at all radii. 

According to these results, dispersion around the predicted scaling rela- 
tion between mass and temperature arises mainly due to the excess of gravi- 
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tational energy present in merging events with respect to the final virialised 
state. This agrees, at least at a qualitative level, with the deviation from 
thermally-supported hydrostatic equilibrium that was seen in Figure 14, IS! 



Observations 



From an observational point of view, several studies have been devoted 



to tes t the mass-temperature relation in clusters of galaxies. IHorner et al 



(1999) claim that the M — T relation is steeper than the traditional scaling, 
following a M oc J 11 - 8-2 - 8 law when the mass is estimated according to the 
/3-m odel. Slopes ste eper than 1.5 are also observed in high-redshift clus- 
ters ( Schindler 1999l .e.g.) and highly-lumino us clusters dEttori and Fabianl 
1999), where isothermality has been assumed. Neumann and Arnaudl (|l999\ 
also using the /3-model to estimate the gravita ting mass, obtain a M — T 



relation consistent with the classical scaling. lAller, et all us.ng 
spatially resolved spectroscopy with Chandra, found consistency with the 
scaling law prediction, but the normalisation was ~ 40% lower with respect 
to numerical simulations. 

However, it should be noted that in the above observational analyses, 
only high-temperature {kT > 4 — 5 keV) clusters are involved. As ha s 



been realized in recent studies (e.g. Mohr et al. 19991 : Ponman et al. 1999), 
the input of energy feedback into the ICM can break the self-similarity 
in the low-tempe r ature clusters and groups. This has been confirmed by 



Nevalainen et al. (j2000h . who inferred a slope a = 1.79 ± 0.14 and a nor- 
malisation significantly lower than the one observed in numerical simula- 
tions, although the classical relation M oc T 3//2 was recovered when the 
lo w-temperature groups are excluded. 

Finoguenov et al. ( 200 ll ) also found that the slope of the M — T relation 
for more massive clusters (i.e. M500 > 5 x 10 13 M Q ) was considerably 



shallower (1.58 ± 0.07) than that obtained for their whole sample (1.78 ± 
0. 09), the normalisatio n being more than 50% lower than the value quoted 
Evrard et al. ( 1996^ . These authors also point out that mass estimates 



m 



based on isothermali ty like the 0- model result in a steeper scaling relation (a 
fact already noted bv lHorner etaDri 999 N l due to the implicit assumption that 
the dark matter density scales as r -2 , thus underestimating the mass at low 
radii and overestimating it at large radii. An additional bias in introduced 
because the value of (3, and hence the asymptotic slope of the gas profile, 
wa s a function of t emperature . 



Xu et alJ (|200lh compared the mass-temperature relation obtained from 



the /3-model mass estimates and an isothermal distribution in hydrostatic 
equilibrium with a NFW potential. They found both models to be in- 
distinguishable, giving a slope steeper than the self-similarity prediction 
unless the clusters at le ss than 3.5 keV are excluded, in agreement with 
Finoguenov et al.1 (121 
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Reference 


A 


M /M vir 







Horner et al. ( 1999) 


200 


0.92 


1 CO 

1.53 


velocity dispersion 




200 


0.63 


1.48 


temperature profiles 




200 


0.51 


1.78 


p-model 




200 


0.47 


2.06 


Sx deprojection 


htton and labian (1999) 


500 


0.40(0.97) 


1.93 


/3-model 


JNeumann and Arnaud (1999) 


1750 


0.75 


1.5 


/3-model 


Mohr et al. (1999) 


1000 


1.15 


1.5 


/3-model 


Nevalamen et al. (2000) 


2000 


0.58(0.75) 


1.77 


temperature profiles 




1500 


0.57(0.73) 


1.77 






1000 


0.52(0.73) 


1.79 






500 


0.41(0.74) 


1.84 




Imoauenov et al. (2001) 


500 


0.65 


1.68 


p-model 




500 


0.77 


1.48 


polyt. /3-model, T > 3 




500 


0.45 


1.87 


idem. T < 4.5 kev 


Xu et al. (2001) 


200 


0.67 


1.60 


/3-model 




200 


0.40 


1.81 


isothermal NFW 


Allen et al. (2001b) 


2500 


0.97 


1.47 


J mw Nl W SCDM 




2500 


0.93 


1.52 


T mw NFW ACDM 


Rttori et al. (2002a N l 


2500 


0.82(1.03) 


1.54 


T mw (NFW or King) 




1000 


0.78(1.18) 


1.76 






2500 


0.88(1.03) 


1.51 






1000 


0.48(1.06) 


1.94 





Table 4.1: Observed M -T relation 



Ettori et al.l d2002al ) have not only studied the slope of the M — T re- 



lation, but they have also extensively investigated the dependency of the 
scaling relations with the dynamical state of the cluster as well as the outer 
radius reached by the observational data. The slope at overdensity 2500 is 
consistent with 1.5, the presence of non-cooling flow clusters significantly 
increasing the scatter. However, the slope at A = 500 steepened up to 
M oc T 12 - 17 ^ - 37 . Regarding the norma lisation, these authors claim that the 
previ ously reported lower values (e.g. iHorner et al.lll999l : iNevalainen et al. 
2000) were due to the steeper slopes found in those studies. 

Observational data are summarised in Table H~T] where the reference and 
overdensity of each study are listed in the first two columns, followed by the 
best-fit values of the normalisation (rescaled to the self-similar prediction) 
and the logarithmic slope a. The number in parentheses corresponds to the 
best-fit Mq when the slope is set to a = 1.5. Particularities of each work, 
such as the assumptions employed to derive the gravitating mass, are noted 
on the right margin of the table. 
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Given the remarkable success of the ACDM model in explaining most 
of the relevant cosmological observations 3 , it is somewhat amazing that ob- 
servers do show suc h a tenacious reluctance to accept it as a fact. Except for 
Allen et alJ (|2001bh . all the mass estimates in Table l4~Tl have been made in 



a $7 = 1, A = 0, h = 0.5 universe. Fortunately, the scaling relation between 
mass and temperature does not depend on cosmology but through the usual 
factor h^ 1 in the mass, which is implicitly taken into account in the units 
we used for M™. 

A more worrying concern is raised by the discrepancy found between the 
masses inferred from different methods and models, as well as that arising 
from different analysis prescriptions. As can be inferred from the compila- 
tion given in Table l4~Tl we are far from reaching an acceptable convergence 
in the observational value neither of the normalisation nor the slope of the 
scaling relation. 



Simulations 

The mass-temperature relation has also been extensively studied by means 
of cosmological numeri cal simulations. I n mo st cases, from the early work 
of Na varro et al. (1995) or Evrard et alJ (|l996h to the recent high-resolution 



simulations of M uanwong et alJ (|2002fi , slopes consistent with the self-similar 
value a = 1.5 are found when only adiabatic physics is considered. 

The normalisation, however, has a little bit more scatter, although not so 
much as that present in the different observational estimates. In both cases, 
this scatter is partly due to the precise value found for the logarithmic 
slope, since the normalisation of the scaling law is closely related to this 
parameter (for a given sample, the higher the best-value of a, the lower the 
corresponding normalisation Mo). 

A compilation of different mass-temperature relations found in numer- 
ical simulations is given in Table 14.21 The structure of the table is the 
same of Table 14.11 in which the observational data were presented. Unless 
otherwise noted, the cumulative mass is compared to the emission-weighted 
temperature Tx- Yoshikawa et al.l rt200(t) gives the normalisation referred to 



both quantities, and lMathiesen and Evrardl ( 200ll ) also computed a spectral 



temperature T s based on a more detailed model of the line emission of the 
ICM gas in the soft X-ray band as well as a subsequent fitting procedure 
closer to the analysis of observational data. 

These authors find that the spectral fit temperature is generally lower 
than the mass or emission-weighted average 4 due to the influence of cooler 



3 As well as som e others... The AC DM model is almost as wonderful as SCDM was ten 
years ago! (see e.g. lCole et"ai]ll994al) 

4 Rather surprisingly, their values of Tx are similar, although slightly lower, than T mw , 
which is reflected in the normalisations quoted in Table H^l This is probably due to the 
presence of cold gas in the most inner part of their clusters. 
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Reference 


A 


M /M vir 


a 




iNavarro et ai. ( iyyo i 


zuu 


1 m 

l.Ul 


1.0 




Hjvraia et ai. (lyyoj 


oc;nn 

zouu 


l.oU 


1.0 






1UUU 


l.OO 


1.0 






OUU 


1 Qfl 

l.oU 


1.0 








1 K 
l.lu 








100 


0.89 


1.5 




Brvan and Norman (1998) 


200 


1.42 


1.5 




Pen CI 9981 


200 


1.19 


1.5 




Eke et al. (19981 


100 


1 


1.5 




Yoshikawa et al. (20001 


100 


1.40 


1.5 


T 

- 1 - mw 






0.91 


1.5 


Tx 


Mathiesen and Evrard (20011 


500 


1.34 


1.52 


T 

- 1 - mw 




200 


1.25 


1.54 






500 


1.85 


1.38 


Tx 




200 


1.65 


1.39 






500 


1.81 


1.48 


T (( 




200 


1.69 


1.51 






500 


1.34 


1.62 


T (1 




200 


1.20 


1.64 




Muanwone et al. ( 20021 


200 


2.28 


1.5 





Table 4.2: M — T relation in previous simulations 
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Relaxed clusters 


Whole sample 




A 


M /M vir 


a 


M /M vir 


a 


J- mw 


zouu 


l.ZU 


l.OD 


l.Uo 


1 71 
1.(1 




i nnn 

1UUU 


1.31 


1.59 


1.30 


1.65 




500 


1.34 


1.53 


1.45 


1.54 




200 


1.29 


1.46 


1.57 


1.33 


T x 


2500 


0.93 


1.61 


0.91 


1.68 




1000 


0.90 


1.62 


1.04 


1.57 




500 


0.84 


1.57 


1.12 


1.40 




200 


0.71 


1.51 


1.15 


1.17 



Table 4.3: M — T relation for our sample 



gas being accreted as part of the hierarchical clustering process. Despite sig- 
nificant departures from isothermality, single-temperature models produce 
acceptable fits for the spectral temperature T s . The unusual coincidence 
that a realistic spectrum has nearly the same shape as an isothermal one 
was explained by these authors in terms of a lack of spectral resolution. 

Although several cosmologies have been assumed in all these works, 
the results do not seem to be very sensitive to this. As pointed out by 
Muanwong et alJ (l2002h . the numerical resolution of the simulation is much 
more important, particularly in those cases where cooling was implemented. 
Increasing the resolution lowers the emission-weighted temperature due to 
the presence of cold, dense gas cores that su rvive in sub-clumps. This has 
been argued by Mathiesen and Evrardl (jioOl) to explain the lower spectral 
temperature detected in minor mergers, proposing the deviations from a 
canonical M — T relation to distinguish these systems from relaxed clusters. 

Table 14.31 summarises our results concerning the mass-temperature re- 
lation, where a separate fit has been performed for relaxed clusters (left 
columns). Our values are consistent with the self-similar prediction in al- 
most every case. Clusters undergoing major merging events are usually 
responsible for the discrepancies. 

The normalisation derived in the present study seems to be in agreement 
with previous numerical work when the mass-weighted temperature is con- 
sidered, but is significantly lower when expressed in terms of the emission- 
weighted temperature. We interpret this result as a resolution effect, maybe 
related to the numerical implementation of gasdynamics (see Section l4.1j) . 
Since Tx is biased towards the highest-density regions of the cluster, and 
both our simulated temperature and density profiles seem to increase con- 
tinuously towards the centre (as expected from the analytical predictions 
based on hydrostatic equilibrium) , the emission- weighted temperature grows 
as smaller scales are resolved in the cluster cores. 
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Figure 4.11: M — T relation for our sample at two different overdensities. Solid lines 
are least-square fits to the relaxed population, dotted lines are used for 
the whole sample. The anal ytical prediction is shown by dashed lines 
and observational data from Alle n et all (JToOlb) by a blue line. 



Indeed, our M — T relation agrees rela tively well with the results o f 
Yoshikawa et al. (|2000h . as well as those of Mathiesen and EvrardI (2001) 



for the mass-weighted temperature. As noted earlier, the values of Tx found 
by these authors are lower than T mw . This is only possible if the densest 
gas is colder than the average. This seems to be the case in their Figure 4, 
where they plot the phase-space diagram of the gas particles of one of their 
clusters. As noted in Section 14.11 we attribute the presence of such cold 
clumps to a numerical artifact associated to SPH softening. 

A possible way to overcome the problem posed by the sensitivity of the 
X-ray related quantities to the detailed structure of the ICM in the central 
regions is to compute the local (instead of average) values at rA- Observa- 
tionally, this has been possible only very recently, thanks to the deprojection 



of high -resolution maps obtained by space-based observatories ([Ettori et al 



2002ah . Although it could be in principle a good way to compare numerical 



and observational results, it is still subject to many uncertainties on the 
observational side. 

Our results at overdensities A = 2500 and A = 200 are compared with 
both observational data and analytical estimates in Figure 14.111 We see 
that the virial prediction based on equation ()4.22|) accurately reproduces 
our simulated M — T relation when only relaxed clusters are taken into 
account. As could be seen in Figure 11.101 merging clusters deviate from the 
expected normalisation at large radii, being more than 50% heavier for a 
given temperature. This biases the best-fit slope towards shallower values 
at r 200) as indicated in Table I4.3( but has an almost negligible effect at 
higher overdensities. 

The best-fit mass-temperature relation found by lAllen et all (|2001bh 
is drawn as a blue line in the left p anel of Figure 14.111 Some authors 
( Ettori et al. 2002al : Mohr et al. 19991 ) find similar normalisations at high 
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over d ensities, whereas others ( Neumann and Arnaudll999l : Nevalainen et al 
2000) quote slightly lower values. The situation is much less clear near the 
virial radius, but there is a general trend in the sense that observational 
estimates seem to indicate hotter temperatures for a given virial mass. 



4.3.2 Luminosity-temperature 



A scaling law that is more directly observable is the correlation between 
bolometric X-ray luminosity and global emission-weighted temperature. It 
has been known long ago that the slope of the relation appears t o fit L x oc T$ 
(e.g. Mitchell et aTlll979l : lEdge and Stewartlll99ll : bavid et alJll99.4 iFabianl 
1994) rather than Lx oc T-^, as would be expected for a simple self-similar 
model, and much theoretical work has been devoted to find the reasons for 
such a discrepancy. 

One of the proposed solution involves some preheating of the intracluster 
gas by a n agent other than gravity, such as su pernova explosions at high 
redshift (jKaiser 199ll : Evrard and Henry 199lh . This additional heating 
would decrease the central density and hence the X-ray luminosity. This 
view was strengthened by the discovery of an appar ent entropy 'floor' in the 



centres of groups and clusters ()Ponman et al.lll999l ). 



However, the amount of heating required is substantial. Although esti- 
mates vary, it seems likely that about 1 keV per particle is needed , which 
may be challenging to explain from supernova heating alone due to an ex- 
cessiv e enrichment of the ICM (|Kravtsov and Yep ell2000l : IValageas and Silk! 
1999). Another difficulty is that observations of the Lya forest indicate a 
much low er temperature for the maj ority of the intergalactic medium at 
z ~ 2 — 3 (jBrvan and Machacekl 2Q00) , a condition that may extend to even 
low redshifts. Although hardly conclusive, these concerns may be pointing 
toward another explanation for these observations. 

Early on, it was suggested that the steepening of the L — T relation 
could be due to s ystematic variation s of the cluster b aryo n fraction with X- 
ray te mperature ( David et al. 1993h . Fabianl ( 1994^ and Markevitch et al. 
(1998) also note that the large scatter in the Lx — Tx plane was mostly 
due to the s trong cooling flow s present in more th an half the low-redshift 



tig 

clusters (e.g. Edge et al.ll992f ). It has been argued ( Allen and Fabianl 1998h 
that these cooling flows, which can account for up to 70% of the total X-ray 
luminosity of a cluster ( Alien! 1998h , could be responsible for the observed 
discrepancy in the slope of this scaling relation. 



Analytical prediction 

Simple arguments based on self-similarity predict that the bolometric X- 
ray luminosity of galaxy clusters should scale with the temperature of the 
gas according to Lx oc T^. Nevertheless, the normalisation of the L — T 
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relation is often treated as a free parameter. In order to attempt a theoretical 
prediction, we will calculate the X-ray luminosity of a spherically symmetric 
cluster of galaxies, assuming that it is dominated by thermal bremsstrahlung 
emission. With this prescription, the cumulative bolometric luminosity can 
be computed as 

L x (r) = r A x [kT(x)} 1/2 Attx 2 dx (4.26) 

Jo V V m P J 



where Ax = 1.2 x 10 24 erg s 1 cm 3 keV" 1 / 2 (|Navarro et alJll99.fl . As 



explained in Section 12.2.^1 this is precisely the definition we use to compute 
the bolometric luminosity of our numerical clusters. 

Although it has already been shown that clusters are definitely not 
isothermal, the term T l l 2 (r) varies very slowly compared to p 2 (r). There- 
fore, at a first approximation we can take its mean value out of the integral 
in (|4.26|) . Since this expression is biased towards the highest-density regions, 
it would be more accurate to use the emission-weighted average. 

After some calculations, we obtain the expression 



Ax rim , ,1/2 Jo'*) 4vrx 2 dx/(f r 3 ) M, 



r 



which can be simplified using the definitions of cumulative baryon fraction, 
overdensity and the structure parameter Q(r) = (p 2 ) / (p g ) 2 . The final 
scaling relation is reduced to the form 



rpA 

1 keV 



(4.28) 

Proceeding in an analogous manner as we did for the M — T relation, we 
express the normalisation Lq in expression ([4.28ft in terms of the quantity 
L o' r = ^ja M o" - 2 - 85 x 1043 h er g s_1 - For the sake of simplicity, we will 
assume a cosmic baryon fraction Q,\ ) jQ, m = 0.14 irrespective of the cosmolog- 
ical model. This is approximately the ratio expected for a ACDM universe 
and it is roughly consistent with the value observed in galaxy clusters (e.g. 
lAllen et a,l.l l2002V assuming a low VLm- 

Another difficulty in comparing observational and theoretical results 
arises from the different scaling with the Hubble constant. Since the ob- 
served luminosity is proportional to cIl(z) 2 , it scales as h~ 2 instead of h. 
For nearby clusters, this effect accounts for a factor of 6 in the normalisa- 
tion depending on whether the comparison is made in units of h = 1, h = 0.5 
or h = 0.7. We take h = 1 as the prescription that is more consistent with 
both observational and numerical methods. 

We test the validity of our analytic estimates of the L — T relation by 
comparing the values of Lx/(f 2 Q) and MT^ 2 . The values of these quanti- 
ties for all our simulated clusters are plotted in Figure 14.121 as well as the 
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Figure 4.12: Normalisation of the L — T relation, corrected for the effects of baryon 
fraction and structure parameter. 



results of equation (|4.28j) for two different over density values. The agree- 
ment is not surprising, since the only assumption we made for deriving the 
analytical prediction was to consider that T 1 / 2 varied much more slowly than 
p 2 ,, which, as we discussed in Section [4.2.21 is a very good approximation. 



Observations 



The luminosity-temperature relation of several cluster catalogs observed 
with the first generation of X-ray satellites showed a departure from the 
predicted slope a = 2 (e.g. lEdge and StewarHll991 : David et al. 1£)93) with 
a measured value of a ~ 3 and a scatter along the mean relation that could 
be reduc ed considera bly once the effect of the cold cores was taken into 
account (|Fabianlll99i ). 

When the impact of cooling flows on the bolometric luminosity and 
emission- weighted temperature of the cluster is properly corrected, the L — T 
relation is shown to be tighter and clos er to (yet not completely consistent) 
the s elf-similar scaling prediction (e.g. lAllen and Fabianl ll 



199. 



Ettori et al 



2001). The scatter is signifi cantly reduced in anal yses that either excise 



the core cooling flow regions ( Markevitch et al. 199ot) or examine samples of 
clusters defined to possess weak cooling cores (jArnanH and Fvrardlll999l) . 



T he steeper L — T dependence found in other studies (e.g. David et al.l 



1993 1: IPonman et al!l99fi : Jcavaliere et al!l997l : IPonman et alll999l : IXue and WiJ 
1200(1 iTozzi and Normal l200ll ) is often justified by invoking some sort of 
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Reference 


A 


L o/ L o 1T 


a 




David et al. (1993) 


200 


0.084 


3.37 




ronman et al. (1996) 


200 


0.468 


8.20 


groups 




200 


0.132 


3.29 


groups+clusters 


All J T7 1 1~ " /inno 1 ! 

Allen and Fabian (1998) 


200 


0.287 


3.22 


coolmg-now (raw) 




200 


1.698 


2.13 


CF (corrected) 




200 


0.191 


2.92 


non-CF 


Markevitch et al. (1998) 


200 


0.414 


2.10 


0.1 - 2.4 keV 




200 


0.347 


2.64 


bolometric 


Arnand and Evrard (1999) 


200 


0.209 


o o o 

2.88 




Xue and Wu (2000) 


200 


0.017 


5.57 


groups 




200 


0.029 


2.79 


clusters 


Allen et al. (2001b) 


o run 

2500 


n one 

0.296 


2.17 


i mw bCDM 




2500 


0.589 


2.08 


T mw ACDM 


Ettori et al. (2002al 


2500 


0.055(0.258) 


2.79 


T 




1000 


0.347(0.677) 


2.37 






2500 


0.076(0.258) 


2.64 






1000 


0.214(0.551) 


2.54 





Table 4.4: Observed L - T relation 



pre-heated of the intracluster gas, which would rise the entropy level in the 
coolest systems. 

A compilation on observational estimates of the Lx — Tx relation is 
given in Table [Ol We see that both the normalisation and the logarithmic 
slope differ considerably from author to author, being extremely sensitive to 
the details of the data reduction and analysis (e.g. instrumental specifics, 
treatment of cooling flows, fitting procedure, etc.). 

Once again, we see that the parameters Lq and a are tightly correlated, 
reflecting the fact that the typical temperature of most cluster samples is 
substantially larger than 1 keV. Thus, the normalisations obtained for a 
self-similar slope a = 2 (quoted between parentheses in Table l4~4"|) are much 
higher than those obtained when the slope of the relation is not constrained. 

Contrary to the M — T relation, in this case the assumed cosmological 
model plays an import ant role in the exa ct value of Lq, while the slope a 
is not so affected (see lAllen et~al~ll2001bl . the only work that included an 
estimate based on a ACDM cosmology). 

The effect of differences in individual gas fractions and structure param- 
eters of clusters o n the observed logarithmic s lope of the L — T relation has 
been examined bv lArnaud and Evrar 3 £299), who came to the conclusion 
that the higher gas concentration in hot clusters is responsible for at least 
half of the steepening of the L — T relation. When cluster mass is estimated 
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Reference A L /L™ a 

Navarro et al. (1995) 200 1.46 2 

Brvan and Norman (4 9981 200 1.27 2.59 

Eke et al. (1998 'l 100 1.08 2 

Bia.lek et a,l. C2001 ) 500 3.85 1.56 core included 

1.25 2.02 core excised 



Table 4.5: L — T relation in previous simulations 



from a /3-model , the temperature dependence of the structure factor Q can 
explain the observed slope. But if virial theorem is used to derive masses, 
then the gas fraction increases with cluster temperature, and it contributes 
as much as Q to the steepening of the logarithmic slope a. 

Simulations 

Despite the significant amount of observational work regarding normalisa- 
tion and slope of the Lx — Tx relation, there are little numerical work done 
in this regard as compared to works done to study the mass-temperature 
relation (some of them listed in Table I4~2*|) . 

This is a reflection of the difficulty in obtaining reliable estimates of the 
bolometric X-ray luminosities. While it is reasonable to assume that the 
results are sensitive to the ICM thermal evolu tion, it is well known (e.g. 
Brvan and Normanlll998l : lYoshikawa et al~ll2000l ) that they are also affected 



by the numerical resolution. 

When equal mass particles are employed in a given simulation, the small- 
est clusters have poorer resolution and hence their luminosities tend to be 
systematically underestimated as compared with more massive ones. 

This effect can steepen the slope of the simulated L — T relation, making 
it closer to the observed one. Although our numerical experiments seem to 
have enough resolution to overcome this problem, it is always necessary to 
keep in mind this consideration when interpreting the results. 

Many of the references cited in the present work did not study the L — T 
relation at all. Among them, only those reporte d in Ta bic 14.5 1 quot e the 
best-fit values found for the parameters Ln and a. Bialek et alJ (j200lh also 
includes a discussion on the analytical derivation of the scaling laws, high- 
lighting the importance of cluster structure, related to Ft, and Q parameters. 

These authors point out that the slope of the power-law fits can be biased 
by a few points at one end. In their case, they had two low-temperature 
systems that appeared to be fortuitously obse rved mergers. In ag reement 



with our results presented in previous chapter iBialek et al.l ()200lh remark 
that these two objects featured steeper density profiles compared to the 
rest of the sample, and used this fact to justify that their luminosity was a 
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Relaxed clusters 


Whole sample 


A 




a 




2500 


1.34(1.11) 


1.89 


0.55 2.86 


1000 


2.35(2.21) 


1.90 


1.24 2.66 


500 


3.44(3.13) 


1.89 


2.23 2.45 


200 


5.57(4.95) 


1.89 


4.25 2.28 



Table 4.6: L — T relation for our sample 



factor of 3 higher than that expected from their gas temperature. Exclusion 
of the core emission, where the core is defined as a circular area of radius 
0.13r2oo> resulted in a best-fit slope consistent with the analytical scaling 
(see Table H31). 

However, this is not what we find in our numerical experiments. In 
Table 14.61 we show the values of Lq and a that best describe the Lx — Tx 
relation of our galaxy cluster catalogue. The best-fit parameters are quoted 
for different overdensity thresholds, making a separate fit for the relaxed 
subset of the sample. The number in parentheses is the normalisation Ln 
when the logarithmic slope is fixed to the canonical value a = 2. 

There is a striking difference between the behaviour of the relaxed clus- 
ters and th e merging sys t ems i n our simulations. Contrarily to the results 
reported bv lBialek et al.l (|200ll ). the latter exhibit significantly lower lumi- 
nosities for a given temperature, mainly due to their lower concentrations 
(and thus, lower central densities). This is not surprising, since, despite 
the steeper inner slope of the density profile, the normalisation p s is signif- 
icantly reduced with respect to relaxed systems (see e.g. the right panel of 
Figure EJJ. 

We also find, in agreement with lBialek et al ] (|200lh . that the low-mass 
end of our sample is dominated by the presence of merging systems, but 
in our case this biases the X-ray luminosity towards lower values, and the 
slope of the L — T relation steepens considerably when low-mass mergers are 
included in a lest-square fit together with the more massive, relaxed clusters 
of galaxies. 

At high overdensities, merging systems raise the slope of the simulated 
L — T relation to values consistent (and sometimes even higher) than those 
observed. However, since these are the least massive of our clusters, caution 
must be exercised because of the lower relative resolution (i.e. less particles 
within the virial radius) in these objects. Nonetheless, a steepening of the 
luminosity-temperature relation would be expected if 



1. Clusters indeed followed a 'universal' baryon fraction profile (as it 
seems to be the case in Figure HI?} . 
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Figure 4.13: L— T relation for our cluster sample. S ymbols are as in previous Figures. 

Observational data from Allen et ail l|200 1bl are shown as solid blue. 
Dotted line is the self-similar relation. Dashed line is the unconstrained 
best fit to relaxed clusters and solid line is the best fit to relaxed clusters 
for a — 2 (see Table |4j 



2. Low-temperature systems were systematically less concentrated. 

In such dependence of the baryon fraction on temperature 

would arise and hence the X-ray luminosity of the low-temperature clusters 
would be significantly depressed with respect to that of the more massive 
objects. This effect would be more evident at high over densities. Near the 
virial radius, the global baryon fraction is expected to be only slightly lower 
than the cosmic value, since gas and dark matter profiles are observed to be 
proportional everywhere except in the core regions. 

The first hypothesis above can be understood as a consequence of hydro- 
static equilibrium and a 'universal' dark matter profile. The second arises 
from the fact that low-mass systems are more prone to merge(in agreement 
with the predictions based on the hierarchical scenario). Therefore, we think 
that this relation between baryon fraction and temperature is not an artifact 
due to a lack of resolution, but has a physical origin. However, this is not 
the only mechanism that can alter the observed Lx — 3x relation. A more 
detailed discussion is given in the next section. 

The steepening of the L—T scaling law can be clearly seen in Figure ll.131 
where the effects of observational 'aperture' (i.e. overdensity) can be appre- 
ciated at first sight by comparing the left and right panels. At A = 2500, 
our results for the whole sample are perfectly compatible with a power law 
as steep as Lx oc T x . 

Observational data from Al len et al. | d2001bh have been represented by 



a blue line in the left panel of Figure 14.131 It is not clear to what extent 
our results are compatible with this fit, but the disagreement is manifest if 
we compare with estimates based on a SCDM cosmology. However, such 
a comparison is not straightforward. These authors, for instance, quote a 
normalisation which is a factor of 2 lower (as well as a steeper slope) for 
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Figure 4.14: Deviations from the self-similar L — T scaling law at different overden- 
sities. Left panel: Baryon fraction. Middle panel: Structure parameter. 
Right panel: Mass-temperature normalisation. 



the same cluster sample if a SCDM cosmological scenario is assumed for the 
analysis. 

Given the many difficulties surrounding the reliable calibration of the 
Lx — Tx scaling law from both theoretical and observational grounds, the 
mass-temperature relation seems to be much more convenient in order to 
use clusters of galaxies as cosmological probes. 



Departure from self-similarity 

We have already shown in Figure R. 121 that our clusters follow quite closely 
the analytical prediction of the L — T relation given by expression (|4.28|) 
when the effects of substructure are taken into account. However, it is 
evident from both Figure I4.1HI and the slopes quoted in Table 14.61 that the 
best-fit Lx — Tx relation deviates significantly from the self-similar scaling 
law L cx T 2 , particularly at high overdensities. 

Since 1)4. 28j) holds, these deviations arise from the factor F^QMq, where 
we recall that denotes the cumulative gas fraction, Q = (p 2 g ) / {Pg) 2 is 
the structure parameter and Mq is the normalisation of the M — T relation. 
It has been often argued that self-similarity implies that these quantities 
must be the same for all clusters, and hence the normalisation Lq would be 
a constant multiple of L™. 

A systematic dependence of Lq on temperature would tilt the observed 
L — T relatio n, changing the best-fit value of the logarithm ic slope a. Several 
authors (e.g. David et al. 199.4 Arnaud and Evrardl 1999h argued that T b 2 ~ 



T would reconcile the slope a ~ 3 derived from observat ional data with the 
self-s imilar prediction. Some attention has been paid (e.g. Arnaud and Evrardl 



1223) 

to the effect of the structure parameter, and almost none to possible 
variations in Mq 
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Figure 14.141 displays the behaviour of these quantities as a function of 
temperature and overdensity. In the top panel, we see that the cumulative 
baryon fraction at r2oo is slightly lower than the cosmic value, and does 
not vary significantly even in clusters undergoing a major merging event. 
However, these systems are much less concentrated than relaxed structures, 
which manifests as a factor of ~ 3 — 4 difference in the structure parame- 
ter, Q. However, this variation is anti-correlated with the normalisation of 
the M — T relation, since Mq tends to increase for lower values of r2oo/r s 
(see Figure H.10|) . These two effects compensate each other and the total 
deviation from the self-similar scaling is relatively small for A = 200. 

As we restrict ourselves to inner radii, both the structure parameter 
and the normalisation of the M — T relation tend to more stable values. 
Now, differences in central concentrations make that smaller systems, whose 
gas distribution is puffer, get lower baryon fractions than relaxed clusters. 
Hence, the bolometric X-ray luminosity is suppressed by a factor of 3 with 
respect to the expected behaviour. As can be seen in Table fPH this effect is 
responsible for a dramatic steepening of the slope of the Lx — Tx relation. 
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Chapter 5 

Stars 



Oh let the sun beat down upon my face 
stars to fill my dream 
I am a traveler of both time and space 
to be where I have been 

- Led Zeppelin : Kashmir (1975) - 



'tat formation and feedback to the interstellar medium (ISM) play 
a key role in the whole process of structure formation. The cosmological 
model provides the space-time framework in which the events take place 
and fixes several key parameters, such as the expansion rate of the universe, 
the amount of matter and the amplitude of primordial density fluctuations. 
The large-scale structure of the universe, the clustering properties of matter 
or the mass distribution in individual objects can be understood in terms of 
hierarchical growth of dark matter fluctuations. Hydrodynamic simulations 
have shown to be a valuable tool to understand the physical conditions of the 
gas at low and mod erate densities. Detailed theoretic al studies of, e.g. the 
Lyman-alpha forest (jCenlll992l : Hernauist et al. 1996) and the intergalactic 



medium could be done thanks to them. 

At higher densities, radiative cooling and star formation become an es- 
sential ingredient in the description of baryons. Unfortunately, our knowl- 
edge of the complex physical processes involved in galaxy formation is much 
less certain. For instance, in the absence of stellar feedback most of the 
baryonic matter woul d have cooled into small mini-g alaxies (|Colel Il99ll : 



Blanchard et al. 1992h . Photoionization and supernova explosions inhibit 



cooling and star formation in very low-mass objects, biasing star forma- 
tion at early times towards high overdensity regions. Feedback has been 
proposed, as we discussed in the previous chapter, to solve the discrepancy 
between the observed and theoretical scaling relatio ns of galax y clusters, as 
well as the problems of angular mo mentum in discs dSilkl 2001 ) and the low 



number of observed dwarf galaxies ( Bullock et al. 2000l ). 



Ill 
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Despite the num erous attempts to include star formation in cosmologi- 
cal simulations (e.g. Springel and Hernquist 2002bl .and references therein), 
progress has been comparatively low with respect to gravity or gasdynamics, 
since the physics involved in star formation and feedback into the interstel- 
lar medium is substantially more complex, and also because it acts on scales 
of a few pc and masses of ~ 10 5 M , which are well below the resolution 
limits. 

The prese nt chapter d e scribe s the results obtained from the model im- 
plemented bv lYepes et al. [see also Section 12. 1.1 j) . The star formation 
history of the universe dAscasibar et al.1 l2002ah is discussed in detail, as 
well as the main differences between the star formation regimes observed in 
clusters and field galaxies. 



5.1 The cosmic star formation history 



The star formation history of the universe, i.e. the global star formation rate 
(SFR) density as a function of redshift, is a crucial test for galaxy formation 
scenarios. Many measurements of this quantity have been undertaken during 
the past few years, and the number and accuracy of available observations 
are still rapidly increasing. 

Although there is a large scatter among different indicators at any given 
redshift, basically all studies find a significant increase (by about one decade) 
in the cosmic SFR densit y from the present d ay to z = 1. However, it is 



still a controversial i ssue (iHopkins etal 
maximum there 



200 If) whether it re aches a broad 

.g.lMarian et alJIl 994 lOisnert et alJbood ) or it declines 



gradually towards higher redshifts (e.g. Sawicki et al.lll997l : Pascarelle et al 
iSlSteidel et~aTlll3. 



From a theoretical point of view, star formation is a highly non-linear 
process, which precludes any simple analytical treatment. Most efforts to- 
wards modeling the star formation history of the universe resort to semi- 
analytical methods or numerical simulations to tackle the feedback mecha- 
nisms that self-regulate the SFR within the hierarchical clustering scenario 
of structure formation. 

Semi-analytical mode ls (see e.g. Kauffmann et al.ll993al : Cole et al.ll994bl : 
Salvador-Sole et al. 199ftl ) construct a large sample of Monte Carlo realiza- 
tions of halo merging histories, using the Press-Schechter formalism. Hy- 
drodynamics, star formation and feedback are implemented through simple 
recipes, whose parameters are fixed in order to match the observed properties 
of real galaxies, specially the lumino sity function and the Tully-Fisher re- 
lation ( Somerville and Primaclj[l9*9 9\ Their high computational efficiency 
allows a fast exploration of the parameter space in search of a viable model, 
but the number of free parameters asks for supplementary modeling based 
on more fundamental physics. 
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Numerical simulations are intended to solve directly the set of differential 
equations of gravitation and gas hydrodynamics, and hence require much 
fewer model assumptions than the semi-analytic approach. However, despite 
the staggering advances in computer technology and numerical algorithms, 
resolution in large-scale hydrodynamical simulations is severely restricted by 
computational and data storage requirements. Although gas can be treated 
self-consistently, some phenomenological recipes are still required in order 
to model star formation and feedback. 

The cosmic SFR densi ty is a direct outcome of semi-an alytical meth- 



ods dSomerville et al. 2001 ) and hydrodynamical simulations dTissera 2000; 



Nagamine et al.112001 ; Ascasibar et al.l2002al : Springel and Hernquistl2 002c^ 



Both of them predict a significant star formation activity at high z, favoring 
the hypothesis of a gradual evolution of the comoving SFR as a function of 
time. 



5.1.1 Observations 

A variety of observations and observational methods have contributed dur- 
ing the past years to our understanding of the evolution of the star forma- 
tion rate density on cosmological scales. Nonetheless, different wavelength 
regimes and estimation prescriptions yield rather discrepant results. In this 
section, the merits and flaws of each spectral range are reviewed, as well as 
the usual procedure to construct a consistent diagram of the SFR density 
as a function of redshift, also known as the 'Madau plot'. 



Wavelength 



Photo metric redshift estimates based on Lyman- limit systems ( St eidel et al 
199qDQ) stimulated the study of the rest frame UV continuum flux around 
2000 A as a preferred wavelength range for objects in the far universe. This 
emission traces the presence of massive (and therefore young) stars and can 
be related directly to the actual star formation rate. It is relatively easy to 
detect at high z as it is redshifted into the optical band. However, older stel- 
lar populations and AGN can make a significant contribution as well, leading 
to an overestimation of the total SFR. On the other hand, dust enshrouding 
star formation regions absorb very efficiently the UV light, re-emitting it 
in the far IR. Unfortunately, there is still some un certainty to account for 
dust extinction and its evolution with redshift (e.g. Adelberger and Steidell 
2000). 



Recombination lines from HII regions can be in principle a more reliable 
estimator of the instantaneous SFR, since the ionizing radiation (A < 912A) 
comes from more massive (younger) stars than the softer UV continuum. 
These optical lines are less dramatically affected by dust extinction and can 
be mapped with higher resolution in the local universe, but they require 
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infrared spectroscopy at moderate redshifts, which has not been possible 



until very recently due to the faintness of the sources involved (|Pettini et al 



1998). For H a , this happens at z ~ 0.5, and for forbidden lines, such as 
[OII]3727A, at z ~ 1.5. There is also a considerable degree of uncertainty 
in the relation between recombination lines and star formation activity. 

Finally, the far infrared (FIR) spectrum of a galaxy can be roughly sepa- 
rated into two components, namely the thermal emission of dust (heated by 
star formation bursts) around A ~ 60 ^m and an infrared cirrus powered by 
the global radiation field, which dominates for A > 100 p, Neglecting cir- 
rus and AGN contamination, the FIR luminosity would also be an excellent 
tracer, but instrumentation is not as developed in the IR and radio bands 
as it is in the optical, and the conversion factor between FIR luminosity and 
SFR is for the moment rather model dependent. 



The Madau plot 

The first step to obtain an observational estimate of the cosmic SFR density 
is the selection of a complete sa mple of galaxies at a given redshift and 
wavelength. Some authors (e.g. Steidel et al. 1999h claim that the HDF 



covers a very small area (~ 5 arcmin 2 ), and hence it is not statistically 
representative of the whole universe. This cosmic variance is difficult to 
quantify, but the strong clustering of high redshift galaxies (and the fact 
that most estimates beyond z = 2 are based on the HDF) seems to indicate 
that this effect could introduce an important systematic bias in the star 
formation rates inferred for galaxies in the high redshift universe. 

Then, the sample must be corrected for incompleteness before construct- 
ing the lumino sity function, which can be done following the V may[ formalism 
(|Schmidtlll968h or re sorting to more elab orate Monte Carlo algorithms (such 
as that proposed by Steidel et al. 1999h . The comoving luminosity density 



is easily obtained by integration over all magnitudes, extrapolating the faint 
end as a Schechter function. The last step is to convert the LF to a star 
formation rate density with the help of a population synthesis model, taking 
into account the absorption of interstellar dust. 

A compilation of recent results is given in Table 15. 1\ where the first 
column indicates the bibliographic reference, as well as the survey from 
which the galaxy sample was extracted. Next columns correspond to the 
chosen wavelength and redshifts of each analysis, and finally the comoving 
luminosity and SFR densities. 

Conversion factors between the last two quantities differ significantly 
from author to author, and in some cases only the luminosity density was 
provided. If there is a SFR estimate in the original paper , then it is quoted 
in Table 15.11 corrected by dust extinction; else, we follow Kennicuttl (1998) 



prescription for an exponential burst and a Salpeter IMF between 0.1 and 
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38.631^8 


-1.97 


-1 


.94 






0.625 


± 





.125 


39.16±g;iS 


-1.44 


-1 


.44 






0.875 


± 





.125 


39.561°;*° 


-1.04 


-1 


.07 


CFRS (Flores et al. (1999)) 


15 fim 


0.35 


± 





.15 


41.97 ±0.25 


-1.46 


-1 


.43 






0.625 


± 





.125 


42.20 ± 0.25 


-1.15 


-1 


.15 






0.875 


± 





.125 


42.53 ±0.25 


-0.82 


-0 


.85 


HDF (Hushes et al. (1998)) 


60 fim 


3.0 


± 


1 


.0 


< 42.26 


-0.98 


-1 


.05 


(Haarsma et al. (2000)) 


1.4 GHz 


0.28 


± 





.12 


26.751°;^ 


-1.17 


-1 


.13 






0.46 


± 





.05 


27.03±g;i* 


-0.89 


-0 


.87 






0.60 


± 





.05 


07 1 o+O- 16 
Z I-XZ_ 26 


-0.80 


-0 


.80 






0.81 


± 





.08 


27.39±g:i| 


-0.53 


-0 


.55 






1.60 


± 





.64 




-0.38 


-0 


.44 



Table 5.1: Observational estimates of the cosmic SFR density at different epochs. /5l 
refers to the luminosity density at the appropriate wavelength, expressed 
in erg s _1 Hz _1 Mpc~ 3 . represents the comoving SFR density (in M Q 
yr _1 Mpc~ 3 ) for our SCDM and ACDM cosmologies. d Original data has 
been corrected for dust extinction. 
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100 M : 

A r M vr- 1 ! -l 1Ax 10_28 L uv [erg s^Hz- 1 ] 

P * [M ° yr ] " \ 1.4 x 10- 41 L 3727 [erg s^z^] ^ 

For dust extinction, we have applied a correction of A(1500-2000A)=1.2 
mag and A(2880A,OII)=0.625 mag, which correspond to factors of 3.02 and 
1.78 respectively. All modifications to the original data are noted by a 
footnote mark on the bibliographic reference in Table I5TT1 



Cosmology 

Most observational values given in Table 15.11 were calculated assuming a 
standard cold dark matter (SCDM) cosmology with dimensionless Hubble 
parameter h = 0.5 and mean matter density O m = 1. Since the conversion 
factor from luminosity to SFR and the correction for dust extinction are 
completely independent of the cosmological scenario, the computation of 
the comoving SFR density in a ACDM (h = 0.7) model only involves the 
transformation of the luminosity densities. 

In a flat universe, the volume enclosed by a solid angle Au> between 
z — Az and z + Az is 

V(z, Az) = ^Au, (d 3 m (z + Az) - d 3 m (z - Az)) , (5.2) 

where d m (z) is the comoving distance. Luminosity scales with luminosity 
distance, d^(z) = (1 + z)d m (z), as d^. Then, assuming that all galaxies are 
located at the center of the interval in z, the luminosity density (and hence 
the total SFR density) will be proportional to 

o ( z ) oc L{Z) oc ^ (5 3) 

MZ) V(z, Az) ^ dl(z + Az) - dl(z - Az) {b - 6) 

The conversion between SFR densities in ACDM and variants of the 
SCDM (h = 0.5) cosmologies can be approximated, to first order in Az, as 



P 



ACDM ^ACDM 



pSCDM /jSCDM 



^n m + n A (i + z )- 3 (5.4) 



Figure ISTl shows this factor, as well as the exact expression for different 
values of Az, which requires a numerical evaluation of the luminosity dis- 
tance in the ACDM model. Although the results are very similar, we have 
chosen to use the exact expression ()5.3|) to compute the SFRs in Table l5~Tl 
in order to minimize errors. 

At low redshifts, the data are consistent for all the independent tracers 
of the star formation activity once the effect of dust obscuration is taken 
into account. The steep increase in the cosmic SFR density from the present 



5.1. THE COSMIC STAR FORMATION HISTORY 



117 




z 

Figure 5.1: Conversion factor between star formation densities in SCDM and ACDM 
cosmological models, for several redshift intervals Az. Dashed line cor- 
responds to the approximation 1)5. 4|) . 

day to z ~ 1 is firmly established. However, there is not yet a full agreement 
on whether there is a characteristic epoch of star formation in the universe 
around z ~ 1.5 or, on the contrary, the SFR declines very slowly towards 
higher z, as more recent observations seem to point out. The results of 
ongoing deep surveys will be extremely helpful in order to reach a definitive 
conclusion from the observational point of view. 

5.1.2 Simulations 

In order to study the cosmic star formation history, several simulations have 
been accomplished using the Eulerian gasdynamical code YK 3 (succintly 
described in Section ^.l.lj) . which incorporates a simple physically motivated 
model to implement star formation and feedback into the ISM. Our aim was 
to investigate the separate effects of cosmology, supernova explosions and 
photoionisation on the self-regulation of the star formation activity. For 
details on the numerical experiments, the reader is referred to Section [7331 

Cosmology 

Three different cosmological scenarios have been considered in the present 
study: the currently-favoured ACDM universe, the 'old-fashioned' standard 
CDM model, and a modified version of the latter: The Broken Scale In- 




Figure 5.2: Evolution of the comoving SFR density in different cosmologies (experi- 
ments ACDM5, SCDM and BSI). The average value over different real- 
izations is plotted as a thick solid line, while the shaded area represents 
the lcr deviation fro m the mean. Dots cor respond to the observational 
data points listed in lAscasibar et alJ (|2002a) 



variance (BSI) mo del which is based on a double inflation scenario (see e.g. 
Kates et al. 199.4 for a description). The post-recombination power spec- 



trum of the BSI model is similar to that of the rCDM scenario, in which 
the dark matter consists of a decaying neutrino. The global star formation 
histories found in our simulations of each one of those three cosmologies are 
compared in Figure IB~2l A crucial factor is the initial normalisation of the 
power spectrum, since this parameter sets the amount and evolution of sub- 
structure, and therefore represents one of the key ingredients in determining 
the rate of cooling and subsequent conversion of gas into stars. 

Cosmological scenarios with CDM alone cannot explain structure for- 
mation on both small and very large scales. However, we have included 
SCDM as a reference model in our comparison. Due to its large power at 
small scales, this scenario produces many structures at high redshift, and 
so it features a flat star formation rate even beyond z ~ 6, with a plateau 
about one decade over the observations (Figure E21 middle panel). It does 
not seem very likely that this difference can be explained by dust absorption 
alone. In this model, too many baryons are transformed into stars, and it 
happens too early with respect to observations. We claim that the cosmic 
star formation history poses an additional problem for the SCDM model. 

On the contrary, the BSI model leads to a clear maximum in the SFR at 
z ~ 2, but the absolute value is too small and the decrease at higher redshifts 
is clearly too steep (Figure f5,2| right). This model has not enough power at 
small scales for producing enough stars at high redshift. As a side effect, the 
cosmic star formation rate density remains too high at low redshifts, since a 
large part of cold gas remains in high density regions triggering further star 
formation activity. 

A cosmological model with non-zero cosmological constant offers the best 
agreement with observations (Figure 15.21 left). Both the overall trend and 
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Figure 5.3: Effects of stellar feedback on the cosmic SFR density, a) ACDM5 sim- 
ulations with supernova feedback parameter A = 200 (solid line) within 
la error (shaded area) and A = (dashed line), b) Same as (a) but for 
ACDM4 simulations run with different values of the overdensity threshold 
T> imposed by photoionization. 

the height of the plateau fit reasonably well the observational data points. 
At present, the AC DM model fits best all known observational data (see e.g. 
Bah call et al 1 ll999h . Our results confirm that this model can also reproduce 
the observed star formation history of the universe when dust extinction is 
taken into account, in which case a significant amount of star formation is 
expected to take place at high redshift. 

A further feature, common to all cosmologies, is the smooth redshift 
dependence of the mean SFR and the small range of scatter. This is in sharp 
contrast to the evolution of individual baryonic clumps or the behavior of 
different regions in a single time-step of the simulation. There, we observe 
huge bursts of star formation (with an increase in activity by a factor of 10 
and higher with respect to the average value), as well as a quiescent regime 
in which the star formation rate remains approximately constant. 

In these simulations, the first stars are created around z ~ 10, which 
seems to be not too unrealistic, but it may depend on resolution. This was 
not a primary issue of the present study. 

Feedback 

Figure shows the results of experiments ACDM5 and ACDM4, from 
where the effects of supernova feedback and photoionization were studied. 
The same trends described in this section were also seen in SCDM and BSI 
simulations. 

Supernovae (left panel of Figure I5.3JI play a key role in the slope of the 
SFR density at low z. Gas heating and evaporation act as a self-regulating 
mechanism that inhibits subsequent gas infall when many new stars are 
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formed. This is specially important in the most massive objects, which show 
a much higher star formation activity when feedback is not implemented. 

In addition, supernova explosions can also blow away the gas reservoirs 
of dwarf galaxies. The effect on the global star formation rate is more notice- 
able at high redshift, when there are still very few galaxies massive enough 
to retain the heated gas within their gravitational potentials. However, 
photoionization (right panel) is the dominant effect at very early epochs, 
preventing the ga s from cooling and forming stars in low density regions 



(|Efstathioul Il992) . As we have discussed in Section I2.1.H the model im- 
plemented in YK 3 takes into account this effect through the overdensity 
threshold parameter T> for star formation. Only in regions denser than T>, 
molecular clouds will be able to form, and the onset of cosmic star formation 
activity is thus retarded until objects become dense enough to screen the 
photoionizing background. 

We found a SFR density consistent with observational data for T> ~ 100. 
This is in roughly agreement w ith results from a more d etailed theoretical 
estimate of the same parameter ((Muecket and Kateslll9 97) , although in that 
case the actual value of T>, as well as the temperature range in which thermal 
instability is an efficient process, depend on the local gas density and the 
intensity of the ionizing flux. 

Resolution and volume 

Star formation and feedback processes are incorporated in gasdynamical sim- 
ulations (either SPH or Eulerian) by means of rather simple recipes which 
pretend to extrapolate the effects of local physics acting on scales several or- 
ders of magnitude below the resolution limit (often at kpc or even higher, see 
Table EUD for these simulations). In the algorithm implemented in YK 3 , part 
of this extrapolation is hidden in two parameters: The supernova feedback 
parameter A and the overdensity threshold for star formation T>. Whether 
this parameterization is a reasonable extrapolation of the underlying sub- 
grid physics depends on the stability of the results against changes in spatial 
resolution of the simulations. 

On the other hand, when one tries to estimate cosmological averaged 
quantities, like the SFR density, it is necessary to check for deviations of 
these quantities due to small number statistics. In the case that matters 
here, a volume averaged quantity, the cosmic variance can have a relatively 
large contribution in the particular determination of the comoving SFR den- 
sity. Thus, we have simulated different volumes with similar resolutions, 
depending on the computational facilities at our disposal. 

To check for these two effects we compare in Figure 15.41 our results for 
ACDM simulations with different resolutions and volumes. The most no- 
ticeable feature that can be observed in this figure is the steeper slope of the 
local SFR for larger volume simulations. As we increase the volume, larger 
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Figure 5.4: Comoving SFR in our 5 ACDM experiments, all with A — 200 and V = 
100. A significant star formation activity is found at high z regardless of 
resolution and volume. This is consistent with the full set of observational 
data points (Table l5"T]l . 



structures are being simulated (i.e. clusters), in which the star formation 
activity is significantly reduced during recent epochs. In simulations with 
smaller volumes, most of the galaxies that form are either isolated or in 
small groups, and hence have a completely different star formation activity. 
The effects of environment on the SFR are investigated in more detail in the 
next section. 

If this trend we found in our simulations is indeed real, we can conclude 
that the effect of cosmic variance might be quite important in the quan- 
titative estimation of the SFR in the nearby universe, and hence it seems 
very likely that it must be taken into account in the determiniation (either 
numerical or observational) of the evolution of the cosmic star formation 
rate density from z > 1 to the present. 

Although we do not find significant differences in our estimates of the 
SFR density from simulations with different spatial resolutions, there is a 
general trend to have more SFR density at early times in the highest reso- 
luti on experiments. T his is a consequence of the overcooling problem (see 
e.g. Balosh et al. 20011 . for a recent revuew), and it is directly related to the 
increasing number of low-mass halos resolved in these simulations, particu- 
larly at high redshift. A detailed analysis of this question would require a 
more realistic treatment of photoionization than the simplistic prescription 
of an overdensity threshold for star formation. 
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Experiment p*(z = 0) p*(z = 3)/p*(0) 



ACDM1 6.50 0.252 

ACDM2 8.60 0.243 

ACDM3 9.51 0.115 

ACDM4 9.01 0.116 

ACDM5 10.9 0.119 

SCDM 12.3 0.115 

BSI 7.37 0.041 



Table 5.2: Present star density (in 10 8 h M Mpc 3 ) and fraction of stars older than 
z = 3 for the different numerical experiments. 

Keeping all these considerations in mind, we should be extremely cau- 
tious when making quantitative estimates of the cosmic SFR density at any 
given z. However, this quantity shows a similar behavior in all ACDM ex- 
periments. This is a very promising result indicating that its qualitative 
evolution (i.e. the shape of the curve) is a robust prediction of our sim- 
ulations. More precisely, we claim that no characteristic epoch of cosmic 
star formation is found for a ACDM universe. The comoving SFR density 
increases with lookback time until it reaches a nearly constant mean plateau 
beyond z = 2. Stars probably formed gradually at high redshift, according 
to most recent observations and taking dust extinction into account. At 
some point around z ~ 1, the global star formation activity would have 
declined sharply until the present day, although individual SFRs of active 
galaxies are kept approximately constant. 

The star density p*(z) shows a similar behavior for both CDM and 
ACDM models. The most remarkable difference is found in the present 
value of /0*(0) (see Table l5~2|) . which can be observationally deduced from 
the local luminosity density, assuming a constant mass-to-light ratio (e.g. 
Marian fit al 111 9961 : iFiikinrita fit aflllflfld ) : 



p° bs (0) ~ 7 x 10 8 h M Mpc" 3 (5.5) 
Based on the fraction of stars in bulges, ellipticals and SOs, Renzinil 



(1998) estimated that about 30% of the stellar content of the universe should 
have already been formed by z = 3, regardless of the cosmological model. 
The amount of stars produced in experiments ACDM1 and 2 (Table l5~2*|) . 
both at z = 3 and z = 0, are in fair agreement with the observationally 
derived values for these quantities. For the other experiments, the fraction 
of stars older than z = 3 is significantly lower, not only beacuse of a low 
star formation activity at high redshift, but also due to a higher SFR in 
the local universe, as well as to the cosmic bias towards low-density regions 
introduced by the small simulated volumes. 
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Th ese results have been recently confirmed by Springel and Hernauistl 
(2002c) who estimated the cosmic star formation history in the ACDM 
model using a similar prescription for star formation and feedback in the 
Gadget code. Through a series of simulations with different resolutions and 
volumes, they claim to have found stable results from the numerical point 
of view. The general features of the cosmic st ar formation history fou nd by 



these authors are similar to those reported in lAscasibar et al. | (|2002ah . Al- 



though sistematically lower, their SFR density is marginally consistent with 
observational estimates. The reason for this discrepacy could be atrributed 
to the extreme feedback model they have used in their simulations. 



5.2 Star formation in clusters of galaxies 



In the hierarchical scenario of structure formation, the first cosmological ob- 
jects are expected to collapse around the highest overdensity peaks, merging 
untily they eventually form the giant ellipticals in the centers of clusters and 
the bulges of the most massive isolated galaxies. These objects contain the 
vast majority of stars older than z = 3, and hence their star formation his- 
tory must be completely different from that of galaxies living in less dense 
environments. 

The strong morphology segregation observed in rich galaxy clusters (e.g. 
lDressler] ll980) testifies the fundamental role played by the environment on 
the evolution of galaxies. Yet, the physical mechanisms responsible for such 
transformations are still a matter of debate. Several processes might alter 
the star formation conditions in cluster galaxies. Some are related to the 
interaction with the intracluster medium (a fact already pointed out by 
Gunn and Gottl ll972) and others account f or the effects o f tidal forces due 
to the gravit ational potential of the cluster jMerrittlll983h or galaxy-galaxy 
interactions ( Moore et al. 199fil . 1998bl . 1999a| ). All these mechanisms can 
produce strong perturbations in the galaxy morphology, which include the 
formation of tidal tails, dynamical disturba nces that appear as asymmetries 
in the rotation curves dDahle et all l200lf> and significa nt removal of the 
diffus e gas ( Giovanelli and HavnesHl985T Valluri and Jogl 1990l : Quilis et al. 
bOQflh 1 . 

Most of these processes are expected to change significantly the star for- 
mation activity of galaxies in clusters. Several studies have addressed the 



1 In this regard, our simulations of clusters of galaxies described in previous chap- 
ters give also an indication that gas stripping is a very important phenomenon to in- 
hibit star formation in cluster galaxies. We do observe that most of the small ha- 
los that orbit the cluster potential do loose their gas very fast. A visual impres- 
sion of this phenomenon can be seen from a series of animations done from one 
of our SPH cluster simulations. They can be accesed from the following web page 
http://pollux.ft.uam.es/gustavo/VIDE0S/LCDM80/clustervideo.html We plan to do 
a detailed quatitative study of this effect in the near future. 
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issue of the influence of the cluster environment on the SFR of disk galaxies, 
but no agreement has been established so far: whereas some authors pro- 
posed similar or even enhanced star formation in cluster spirals than in the 
field, some others claim quenched SFRs in cluster spirals. This discrepancy 
could arise from non-uniformity of the adopted methods and wavelengths or 
from real differences in the studied clusters. 

If the star formation rates are reduced when galaxies are accreted into 
larger groups or clusters, the cosmological implications could be profound, 
since approxima tely 30 — 40% of the p resent-day galaxies are observed not to 
be isolated (e.g. ICarlberg et al.l l2001.and refere nces therein), in agree ment 
with the results found in numerical models (e.g. iGottlober et ajJl200lh . As 
structure builds up in the universe, more and more galaxies can be found in 
groups and, if these environments serve to terminate star formation, it will 
be clearly reflected in the evolution of the universe as a whole, explaining 
part of the observed decline in the global SFR at recent epochs. 



5.2.1 Observations 



The effect of local environment on galaxy evolution in general is not well 
understood. Many observational studies of environmental effects have been 



dramatically from more common galaxies (e.g. Dressier 1980: Dressier et al. 


198,4 


Couch and Sharpies 198' 


7; Baloeh et, al.lll997l. Il99fll: IPoggianti et al. 


x y tJtJ • 


Moss and Whittk 


2000" ICouch et al.l 12001: Solanes et all 12001). In 


particular, balogh et al. 


(1997, 


1998) showed that the mean cluster galaxy 



star formation may be supressed as far as twice the virial radius from the 
cluster center, relative to a field sample with similar bulge-to-disk ratios, 
physical disk size, and luminosity, which suggests that the observed decrease 
in the SFR may not be fully explained by the morphology-density rel ation. 
This conclusion has also been suggested by ( Hashimoto et al.lE9 98) , who 



studied the relation between local galaxy density and star formation activ- 
ity in the Las Campanas Redshift Survey, and similar results have been 
obtained in more re cent studies based on the 2dF Galaxy Redshift Survey 
( Lewis et al.1l2002bl) and th e Early Data Release of the Sloan Digital Sky 



Survey ((Gomez et al.l l2002). 



From a wide field ph otometric survey using the SUBARU telescope, 
Kodama and Smai 3 (|200ll ) reported the detection of an abrupt change in 
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the colours of galaxies at a 'critical' projected surface density of ~ 200 h 
gal. Mpc~ 2 . This corresponds to a radial distace of ~ 1 Mpc from the 
cluster ce nter. Such a 'break' has been confirmed by the SFRs derived from 



the SDSS lGomez et al 



(2002), although it occurs at an order of magnitude 
lower surface density, i.e. at a larger clustercentric radius. Further observa- 
tions are needed in order to discern whether this difference is physical (i.e. 
due to the differences in redshift or luminosity limits of both studies) or an 
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artifact due to the specifics of the analysis techniques. 

In any case, the 'critical' density (or radius) where the star formation 
rate of galaxies is found to change from that of the field seems to be well 
beyond the virialised regions of clusters of galaxies. Since the mechanisms 
that alter the SFR in the cores of rich cl usters, e.g. ram-press ure stripping 
dGunn and Gottl 19721 : Quilis et"al~l 200(1). galaxy ha rassment, ( Moore et al 



1999a]), or tidal disruption (|Bvrd and Valtonenlll99ol '> . are unlikely to be im 



r 



portant at low overdensities, these results point in the direction that the SFR 
is significantly reduced due t o close interactions in less massive groups within 
the infall region of clust ers ( Zabludoff et al. 199fil : Zabludoff and Mulchacv 



1998; 



Gomez et al 



2002l .e.t: 



In this scenario, environmental effects associated to hierarchical structure 
formation could be responsible for some of the observed decline in the cosmic 
SFR from z ~ 1 to the present. Results from the SDSS jGomez et al.ll2002l ) 
seem to support this view, indicating that the break in the SFR-density 
relation corresponds approximately to the turn-around radius of the studied 
clusters (a few times ?*2oo)) which roughly marks the limit of the gravitational 
influence of the dark matter halo, but is too large for the ICM gas to have 
any significa nt influence on the in falling galaxies. 

Although Eewis et all (2002b) found that the relation between SFR and 
local density was independent on the cluster velocity dispersion and presum- 
ab ly mass, a spectral ana lysis of galaxy groups in the 2dFGRS accomplished 
bv lMartinez et al. (2002) concludes that even the environment of low mass- 
systems (M ~ 10 13 Mq) is effetive in diminishing the process of star forma- 
tion in their member galaxies, but the fraction of spectral types associated 
to star forming gal axies decreases with group virial mass. Furthermore, 
Balogh et al. (|2002h have also found some evidence that the morphology- 



density relation is correlated with the cluster X-ray luminosity. 

To evaluate the co ntribution of clusters o f galax ies to the current star 
formation rate density, llglesias-Paramo et HI {2002) computed the Ha lu- 
minosity function of the central regions of two nearby clusters, obtaining 
the total SFR per unit volume of the member gala xies. Multiplying by the 
cluster density in the universe ( Bramel et al. 2000h . these authors conclude 
that the SFR in type 2 and type 1 clusters can account for approximately 
0.25% and 10.8% of the cosmic stellar production at z = 0, respectively. 



5.2.2 Simulations 



The environmental effects on the star formation activity that can be seen 
in numerical experiments are due to differences in merging history and re- 
moval of the hot gas r eservoir that surrounds every isolated galaxy. In semi- 
analytical model s (e.g. Kauffmann et al1ll993bl : Somerville and Primackfl 999 : 
Gole et alJEoonh . it is assumed that this supply of fresh fuel for star forma- 



tion is lost when galaxies are accreted into larger haloes. Therefore, in this 
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simple scheme, star formation rates begin to decline for any satellite galaxy, 
whether in a poor group or in a rich cluster. 

Although these models are able to successfully match some of the ob- 
served tr ends, such as the radial gradients in SFR within the virial radius of 
clust ers (|Balogh et al.ll20nol : Ipiaferio et al.ll200ll : Oka moto and Nagas himal 
they l ack at present the capability to explore how galaxies lose their 



Bekki et al.l EoQ2 N ). One strength of numerical simulations is 



gas (see e.g. 

that they address the problem from first principles, making it possible to 
constrain the contribution of the proposed mechanisms acting on the SFR, 
as well as the role played by different environments in galaxy formation and 
evolution. 



Oottlober et al.l ((,2001) have studied the merging history of dark matter 



halos as a function of their environment, finding that haloes located inside 
clusters have formed earlier than isolated objects of the same mass. More- 
over, they showed that at higher redshifts (z ~ 1 — 4), progenitors of cluster 
and group galaxies have 3-5 times higher merger rates than isolated galax- 
ies. Mergers of galaxies are thought to play a crucial role in the evolution 
of the SFR. In particular, the inflow of material may serve as a source of 
high-density gas and therefore increase the star formation activity. 

Amongst all the numerical experiments we have accomplished with YK 3 , 
only ACDM 1 and 2 sample a volume large enough to contain clusters and 
groups of galaxies. 

In Figure 1531 we plot the contribution to the cosmic SFR density due to 
the progenitors of present-day objects, according to the mass at z = 0. The 
stellar production of objects that end up in haloes more massive than 10 13 
M© at present (clusters, large groups) are plotted as solid lines, whereas 
the SFR density of those that are in less massive objects at the end of the 
simulation (isolated galaxies or small groups) is shown by the dashed lines. 
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The most massive objects have their peaks of star formation activity 
around z ~ 3, undergoing passive evolution since z ~ 2. This is in excellent 
agreement with the merging history found bv lGottlober et al.l (|200lh from N- 
body simulations. The progenitors of these objects were responsible for the 
bulk of star formation at high redshift, and could be very likely associated 
with the Lyman Break Galaxy population. 

On the other hand, galaxies that are isolated at z = or galaxies residing 
in small groups tend to be m uch younger. In agreement with the results of 
llglesias-Paramo et al. these objects are responsible for most of the 



current stellar production. In these low-density areas, halos massive enough 
to retain their gas content display a nearly constant SFR with sporadic burst 
episodes. 

Unfortunately, the simulated volume in these numerical experiments is 
too small to be cosmologically representative. In fact, the larger volume 
covered by ACDM1 contains only one galaxy cluster of approximately 3 x 
10 14 Mq, whereas the simulation ACDM2 contains four massive groups of 
about 5 x 10 13 Mq. Comparing the solid lines in the left and right panels 
of Figure 15.51 we find qualitatively the expected scenario: the peak of star 
formation in the progenitors of the more massive cluster happened earlier 
(at z ~ 4.5) than for the less massive groups (z ~ 3). This lends some 
support to the claim that the effect of environment on the S FR of infalling 
galaxies is related to the group mass ( Martmez et al.ll2002l ) or luminosity 
(|Balogh et alll2002h . 



As we have seen in Section I5.1.2| photoionization and supernova explo- 
sions inhibit cooling and star formation in very low-mass objects, provid- 
ing a physical mechanism to suppress the over-cooling of gas at very early 
epochs. Thus, the formation of the first stars is biased towards high over- 
density regions that will collapse later to form the cores of the present-day 
galaxy clusters. As their hierarchical assembly proceeds, the gas component 
of these objects is heated by shock waves produced during mergers, and 
the cooling times become too long for efficient star formation to take place. 
Consequently, the SFR drops drastically once the initial reservoir of cold 
gas has been exhausted. 

At the same time, new galaxies will form and produce stars. Some 
of them will remain in rather isolated areas until present (dashed lines in 
Figure 15.5(1 and some others will fall recently into the potential wells of 
the clusters. It seems very likely that tidal inte ractions trigger intense s tar 
formation bursts in these infalling galaxies (e.g. Laverv and Henrv 1$88). 

Observations and simulations coincide in pointing out that star for- 
mation activity in galaxies located in dense cluster environments decays 
(or even stops). T h us, th eir colours become redder, giving rise to the 



Butcher and Oemlerl (|1978l ) effect. Moreover, these objects are expected to 



experience many close encounters and sometimes mergers with neighpour- 
ing cluster galaxies, which could explain their morphological transformation 
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from an initial spiral Hubble type into t he SO population tha t is observed 
to dominate the inner regions of clusters (Dr essier et al .11 199 71 ). In fact, our 
high resolution simulations of cluster formation done with GADGET show 
that dark matter sub-haloes orbit many times in the cluster potential, loos- 
ing most of their gas envelope, but with keeping most of their dark matter 
content 2 . 



2 see http://pollux.ft.uam.es/gustavo/VIDEQS/LCDM80/clustervideo.html for de- 
tailed animation of dark subhalo orbiting within the cluster potential 



Chapter 6 

Conclusions and future 
prospects 

So here we are again to experience the bitter, scalding end 
and we 're the only ones who can perceive it 
But others sing of beauty and the story that 's unfolded 
as one that deserves praise and ritual 

- Bad Religion : Pessimistic Lines (1988) - 



6.1 Conclusions 

/ ^PtOtll a set of numerical experiments, using different codes and initial con- 
ditions, we have been able to compile a database of numerical galaxy clusters. 
High resolution (spatial and mass) has been achieved in selected areas of the 
computational box by means of refinement techniques. The simulations were 
carried out within the framework of the ACDM cosmological scenario, which 
is the most favoured by current observational data. The initial realisation 
was equivalent to a set-up with 1024 3 homogeneous particles. 

From our cluster database, we have been able to study the physical 
properties of both the dark and baryonic components. Most of the work 
accomplished for the present thesis focuses on the analysis of the simulated 
clusters at z = 0. Below we summarise the major conclusions that can be 
derived from the results of this analysis. 

6.1.1 Dark matter 

1. Numerical predictions about the dark matter distribution are robust 
within the resolution limits. ART and Gadget give consistent results. 
Perturbations to the dark matter potential due to the presence of 
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baryons are negligible on cluster scales when only adiabatic physics is 
taken into account. 

From a sample of 15 clusters, 6 have been catalogued as relaxed sys- 
tems, 4 as minor mergers and 5 as major mergers. The amount of 
dynamical activity is significantly higher in the low-mass end of the 
sample, in agreement with the expectation from hierarchical assembly. 

There are systematic differences between the inner profiles of our nu- 
merical clusters according to their dynamical s tate. Relaxed c luster s 
are better described by the profile proposed by Navarro et alJ dl997l'l. 



while m inor mergers are closer to the form advocated bv lMoore et al 



( 1999bh . Vlajor mergers follow a pure power-law distribution for more 



than one order of magnitude in radius; their profiles are even steeper 
than p oc r~ 15 . 

4. The logarithmic slope of the mass profile of relaxed haloes does not 
reach any asymptotic slope at our resolution limit. Indeed, it seems to 
increase above the NFW value for small enough radii. More resolution 
is needed in order to verify this conclusion. 

5. Merging systems are systematically more extended than relaxed clus- 
ters. This effect completely disrupts the usual c — M200 relation, al- 
though some correlation is still found in the p s — r s plane. The best-fit 
values of these parameters are consistent with those inferred from ob- 
servations. 

6. The phase-space density of our clusters is well fitted a power- law 
p/a 3 = 2.24 x 10 4 r -L875 , similar to that found by Tavlor and Navarrcl 
(2001) for galaxy-size haloes. 



7. We find no clear evidence of a 'universal' profile of the anisotropy 
parameter (3=1 — ai/o^.. Although there are hints of some average 
trend, the scatter shown by individual profiles is extremely high. An 
isotropic velocity field (j3 = 0) can be confidently ruled out in relaxed 
clusters and minor mergers. 

8. Angular momentum grows roughly proportional to r is our clusters, 
in agreement with previous numerical work. Rotational support from 
bulk rotation is negligible with respect to the contribution of random 
motions. The eccentricity of individual particle orbits is approximately 
e = 0.5 throughout the whole cluster, but two of our objects deviate 
from this behaviour. 

9. The inclusion of angular momentum and the choice of the initial condi- 
tions are the most delicate issues concerning the spherical infall model. 
Our clusters are very well described by ~ 3a peaks of the density field, 
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smoothed on ~ 1 h Mpc scales, when the eccentricity of the orbits 
is fixed to e = 0.5. 



6.1.2 Gas 



We find additional evidence of non-physical entropy losses in the stan- 
dard implementation of the SPH algorithm. This effect can be of crit- 
ical importance in the computation of the entropy profiles in groups 
and clusters of galaxies . The new formulation of SPH proposed by 
Springel and seems to overcome this problem; high- 



resolution simulations accomplished with this code are consistent with 
the results of the Eulerian code ART. An excessively low number of 
particles can produce an isentropic core due to numerical effects. 

2. The assumption of hydrostatic equilibrium is justified for the clusters 
classified as 'relaxed' on dynamical grounds, although some deviations 
occur in the outer parts of the cluster. It holds only marginally for 
minor mergers, and not at all for systems undergoing a major merging 
event. 

3. A polytropic equation of state with 7 ~ 1.18 provides a good descrip- 
tion of the ICM gas. Only major mergers deviate from a polytropic 
relation. 

4. Under the assumptions of hydrostatic equilibrium, polytropic equation 
of state, a universal profile for the dark matter, and constant baryon 
fraction at large radii, the temperature profile is determined by the 
underlying dark mat t er dis tribution. Our results are consistent with 
Komatsu and Seljak (|200ll ). who follow a slightly different prescrip- 



tion. The temperature profile derived from a NF W profile provides an 
excellent fit to the simulation data. 

5. Projected X-ray temperature profiles are very similar to those observed 
in real clusters. Although the situation is less clear in the inner regions, 
probably in connection with non-adiabatic phenomena, an isothermal 
temperature profile can be ruled out both on observational and theo- 
ret ical grounds . Our data is very well fitted by the formula proposed 
by Loken et al. (|2002h . but no agreement is found regarding the best- 



fit values of the free parameters. 

The density profile of the ICM gas can also be derived by the pro- 
cedure outlined above (4). Again, it is entirely determined once the 
parameters describing the dark matter halo are known. Nevertheless, 
the agreement with numerical data is considerably poorer than for the 
temperature profile. 
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7. The relationship between dark and baryonic profiles is reflected also in 
the scaling relations between global quantities. The M — Tx relation is 
predicted to follow the self-similar scaling relation, particularly at high 
over densities. As we go to larger radii, the density of merging systems 
is significantly higher. This adds to the fact that our prescription 
severely underestimates the gas density at large radii. 

8. Our mass-temperature relation agrees relatively well with previous in- 
dependent numerical experiments. However, some discrepancies are 
found in the emission-weighted temperature due to differences in the 
inner shape of the profile. We claim that this is probably a resolu- 
tion effect. The agreement with observational estimates based on the 
ACDM model and a NFW profile is good; other observational methods 
tend to predict slightly lower normalisations, but the scatter is rather 
large. 

9. The Lx — Tx relation can be derived from the definition of X-ray 
luminosity, substituting the appropriate profiles for the gas density and 
temperature. The scaling law so obtained is not expected to follow the 
self-similar prediction Lx oc T x , since it depends on three additional 
factors: the square of the baryon fraction, F 2 , the structure parameter, 
Q, and the exact normalisation of the M — T relation, Mq. When these 
factors are taken into account, the agreement between predicted and 
simulated scaling laws is extremely good. 

10. At low overdensities (e.g. A = 200), the baryon fraction is close to 
the cosmic value for all clusters, the structure parameter Q increases 
dramatically with Tx (i.e. low-mass haloes are much less concentrated 
due to merging), but Mq is a decreasing function of temperature. The 
last two effects roughly compensate, so they simply introduce some 
scatter in the Lx — Tx relation observed at the virial radius. 

11. At high overdensities (e.g. A = 2500), both the structure factor and 
the normalisation of the M — T relation are almost the same for every 
cluster, but the baryon fraction is lower in the less-concentrated low- 
temperature clusters. Hence, the X-ray luminosity of these systems 
is suppressed with respect to that of the massive clusters, and the 
observed Lx — Tx relation steepens considerably. 

12. Our normalisation of the Lx — Tx relation is significantly higher than 
that found in previous numerical studies. This could well be a nu- 
merical artifact due to our higher resolution. Low-resolution studies 
would not be able to solve the high-density regions near the centre of 
the cluster, and hence the bolometric luminosity would be underesti- 
mated. This issue must be investigated in more detail. 
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13. Observational estimates also yield lower luminosities for a given tem- 
perature. The discrepancy amounts to a factor of 2 for the only esti- 
mate based on a ACDM cosmology, and it is much worse for SCDM 
estimates. However, the comparison between different cosmologies is 
not straightforward. 



6.1.3 Stars 

1. A ACDM cosmology shows the best agreement with observations of 
the cosmic star formation history. Standard CDM tends to overpredict 
the cosmic star formation rate density, while BSI fails to form enough 
stars at high z due to its reduced power at small scales. 

2. Photoionisation fixes the onset of star formation in the universe, sup- 
pressing star formation in low-mass objects at high z. Supernovae 
feedback is very efficient to self-regulate the process of converting gas 
into stars, preventing the most massive halos to form too many stars 
at later times. 



3. In contrast to the first observational estimates dMadau et al.l[l9 96V 



but in agreement with more recent data ( Steidel et alJ 199£ ) , the 



ACDM model predicts almost no drop in the cosmic SFR for 2 < z < 5. 
Star formation seems to be a gradual process with no characteristic 
epoch. 

4. The comoving simulated volumes in our experiments are comparable to 
those covered by observations. Thus, a non-negligible error is expected 
in the observational measurements due to cosmic variance associated 
with small volume statistics. This is most important when one wants 
to determine the steep slope of the SFR density from z ~ 1 to the 
present time. 

5. Star formation histories are very different for objects that end up in 
the cores of clusters and isolated galaxies. Galaxies in cluster cores 
have much older stellar populations, and their activity is drastically 
reduced when they loose their cold gas reservoirs. 

6. In agreement with recent observations, we find that the mechanism 
quenching the star formation in infalling galaxies is efficient in low- 
mass groups (M ~ 10 13 Mq) as well as in clusters. 



6.2 Future prospects 

Our database still has a lot of juice to be squeezed from it. In particular, 
we have restricted our analysis to the present epoch, and there is a lot of 
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work that can be done regarding the study of the evolutionary trends of 
both dark matter and gas properties. 

These numerical simulations are also very interesting to investigate re- 
alistic cluster mergers. All our objects have experienced a major merging 
event since z=0.6. With an appropriate time resolution, we will be able to 
study in detail the physics of mergers and put constrains in the relaxation 
time-scale as a function of the merger characteristics. Moreover, we can 
simultaneously compare the X-ray properties of the system with the dark 
matter distribution before, during and after the merger. Then, an unbiased 
tracer of the dynamical state of a cluster could be applied to real systems in 
order to understand the combined observations in X-ray and optical wave- 
lengths. 

On the other hand, new high-resolution simulations are being carried 
out at the present time. Their aim is to add more clusters to the present 
sample, not only to improve statistics, but also to enlarge the mass coverage. 
This is particularly important for the study of the scaling relations found 
between dark matter and gas properties, as well as to test the validity of 
our prescription relating the structure of both components. 

Finally, the implementation of cooling and star formation would offer 
new insights on the process of galaxy formation and evolution. Some steps 
are also being taken in this direction. In the long term, we would like to 
self-consistently model photoionisation and chemical enrichment, computing 
the intensity of the ionising background at each timestep, as well as metal 
production and advection. 

From the point of view of analytical modelling, we would like to extend 
our treatment of spherical collapse in order to predict the evolution of the 
dark matter profiles. It would be also extremely interesting to study the 
multiplicity function of dark matter haloes in combination with the statistics 
of Gaussian random peaks. Regarding the gas, taking into account infall 
and rotational motions into the hydrostatic equilibrium equation will surely 
improve our theoretical estimations of the ICM density and temperature 
profiles, as well as the normalisation of the scaling relations. 
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6.3 Conclusiones y perspectivas futuras 

6.3.1 Materia oscura 

1. Las predicciones numericas sobre la distribution de materia oscura son 
robustas dentro de los Hmites de resolution. ART y GADGET dan 
resultados consistentes. Las perturbaciones del potential de materia 
oscura debido a la presencia de bariones son despreciables a la escala 
de los cumulos, al menos cuando solo se considera la fi'sica adiabatica. 

2. De una muestra de 15 cumulos, 6 han sido catalogados como sistemas 
relajados, 4 como fusiones' menores y 5 como fusiones mayores. La 
magnitud de la actividad dinamica es significativamente mayor en el 
extremo de baja masa de la muestra, como se podia esperar debido a 
la organization jerarquica. 

3. Hay diferencias sistematicas entre los perflles internos de nuestros 
cumulos numericos dependiendo de su estado dinamico. Los cumulos 
relajados se describen mejor con el perfll propuesto por Navarro et 
al. (1997), mientras que las fusiones menores estan mas proximos a la 
forma defendida por Moore et al. (1999). Los fusiones mayores siguen 
una ley de potencias pura para mas de un orden de magnitud en el 
radio; sus perflles tienen una pendiente incluso mayor que p oc r~ 15 . 

4. La pendiente logaritmica del perfil de masas de los halos relajados no 
alcanza ninguan pendiente asintotica para nuestro lmiite de resolution. 
De hecho, parece aumentar por encima del valor de "NFW" para radios 
sufitientemente pequenos. Se necesita mayor resolution para verificar 
esta conclusion. 

5. La extension de los sistemas fusionados es sistematicamente mayor que 
la de los cumulos relajados. Este efecto distorsiona completamente la 
relation c — M200 usual, aunque todavia se puede encontrar algo de 
correlation en el piano p s — r s . Los valores que mejor ajustan estos 
parametros son consistentes con los inferidos a partir de las observa- 
tions . 

6. La densidad del espacio de fase de nuestros cumulos se ajusta bien con 
la ley de potencias p/cr 3 = 2.24 x 10r~ L875 , similar a la encontrada 
por Taylor y Navarro (2001) para halos del tamano de galaxias. 

7. No encontramos una evidencia clara de un perfil universal en el para- 
metro de anisotropi'a, (5 = 1— < Oq/o^. Aunque hay indicaciones de 
una tendencia promedio, la variabilidad de los perflles individuales es 
extremadamente alta. Podemos descartar con seguridad la existencia 
de un campo de velocidades isotropo ((3 = 0) en los cumulos relajados 
y en los fusiones menores. 
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8. El momento angular crece aproximadamente de forma proporcional a 
r, consistentemente con trabajo numerico previo.El soporte rotacional 
es despreciable con respecto a las contribuciones debidas al movimiento 
aleatorio. La excentricidad de las orbitas de partfculas individuales es 
aproximadamente e = 0.5 para todo el ciimulo. Dos cumulos se desvian 
de este comport amiento. 

9. La inclusion del momento angular y la election de las condiciones ini- 
ciales son los elementos mas delicados en lo que respecta al modelo 
de caida esferico. Nuestros cumulos se describen muy bien con picos 
w 3a del campo de densidades, suavizados en las escalas f» l/i -1 Mpc, 
cuando se fija la excentricidad de las orbitas a e = 0.5. 

6.3.2 Gas 

1. Encontramos evidencia adicional de perdidas de entropfa no ffsicas 
en la implementation estandar del algoritmo SPH. Este efecto puede 
tener una importancia critica en el calculo del perfil de entropfa en 
grupos y cumulos de galaxias. La nueva formulation de SPH propuesta 
parece solucionar este problema; Las simulaciones de alta resolution 
realizadas con este codigo son consistentes con los resultados del codigo 
Euleriano ART. Un mimero excesivamente bajo de partfculas puede 
producir un micleo "isentropico" debido a efectos numericos. 

2. La hipotesis de equilibrio hidrostatico esta justificada para los cumulos 
relajados desde un punto de vista dinamico, aunque se tienen algunas 
desviaciones en las partes externas del ciimulo. Es valida solo de forma 
marginal para fusiones menores y falla completamente para sistemas 
en proceso claro de fusion. 

3. La ecuacion de estado del gas intracumular se puede describir bien 
mediante un politropo de mdice 7 ps 1.18. Solo las fusiones mayores 
se desvian de una relation politropica. 

4. Bajo las hipotesis de equilibrio hidrostatico, una ecuacion de estado 
politropica, un perfil universal para la materia oscura y una fracion 
barionica constante a radios grandes, el perfil de temperatura queda 
determinado por la distribution de materia oscura subyacente. Nue- 
stros resultados son consistentes con los de Komatsu y Seljak (2001), 
que tienen un modelo ligeramente diferente. El perfil de temperatura 
obtenido a partir de un perfil de NFW proporciona un ajuste excelente 
de los resultados de la simulacion. 

5. Los perfiles de temperatura de rayos X proyectados son muy similares a 
los observados en cumulos reales. Aunque la situacion es menos clara 
en las regiones internas, lo que esta probablemente relacionado con 
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fenomenos no adiabaticos, se puede descartar un perfil de temperaturas 
isotermico por razones tanto teoricas como observacionales. Nuestros 
datos se ajustan muy bien con la formula propuesta por Loken et al. 
(2002), pero no se encuentra concordancia segiin los valores de mejor 
ajuste de los parametros libres. 

6. El perfil de densidad del gas ICM se puede obtener tambien con el 
procedimiento descrito mas arriba. De nuevo, queda completamente 
determinado una vez que se conocen los parametros que describen el 
halo de materia oscura. De todas formas, la concordancia con los datos 
numericos es considerablemente peor que con el perfil de temperatura. 

7. La relation entre los perfiles oscuro y barionico se refleja tambien en 
las relaciones de escala entre los promedios globales. Se predice que la 
relation M—T x sigue una relation de escala de autosimilaridad, en par- 
ticular, para sobredensidades altas. Cuando incrementamos el radio, 
la densidad de los sistemas que se fusionan se incrementa significati- 
vamente por encima del valor esperado por su temperatura. Esto se 
acumula con el hecho de que nuestro modelo infraestima severamente 
la densidad del gas para radios grandes. 

8. Nuestra relation temperatura-masa concuerda relativamente bien con 
experimentos numericos previos independientes. Sin embargo, se en- 
cuentran algunas discrepancias en la temperatura pesada con las emi- 
siones debido a diferencias en la forma interna del perfil. Pensamos 
que es probablemente un efecto de la resolution. La concordancia 
con estimaciones observacionales basadas en el modelo ACDM y el 
perfil NFW es buena; otros metodos observacionales tienden a prede- 
cir normalizaciones ligeramente menores, aunque con una variabilidad 
significativa. 

9. La relation L x — T x puede obtenerse de la definition de la luminosidad 
de rayos X, substituyendo los perfiles adecuados de densidad de gas y 
temperatura. No se espera que la ley de potencias que se obtiene asi 
siga la prediction autosimilar L x aT x 2 , porque depende de tres factores 
adicionales: el cuadrado de la fraction barionica, F 2 , el parametro de 
estructura, Q, y la normalization exacta de la relation M — T, Mq. 
Cuando estos factores se toman en cuenta, la concordancia entre las 
leyes de escala predecidas y simuladas es extremadamente buena. 

10. Para sobredensidades bajas (e.j. A = 200), la fraction barionica es 
cercana al valor cosmico para todos los cumulos, el parametro Q de 
estructura crece dramaticamente con Tx (es decir, los halos poco ma- 
sivos estan mucho menos concentrados debido a la fusion), pero Mq es 
una funcion decreciente de la temperatura. Los dos ultimos efectos se 
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compensan aproximadamente, asi que introducen alguna variabilidad 
en la relacion L x — T x observada en el radio virial. 

11. A grandes sobredensidades (e.j. A = 2500) el factor de estructura y 
la normalizacion de la relacion M — T son la misma para todos los 
cumulos pero la fracion barionica es mas baja en los cumulos menos 
concentrados y con baja temperatura. Por lo tanto, la luminosidad 
en rayos X de estos sistemas es menor respecto de la de los cumulos 
masivos, y la relacion L x — T x observada aumenta su pendiente con- 
siderablemente. 

12. Nuestra normalizacion de la relacion L x — T x es significativamente 
mayor que la encontrada en estudios numericos previos. Esto se podrfa 
deber a un artefacto numerico de nuestra mayor resolucion. Los es- 
tudios a baja resolucion podrian ser incapaces de resolver las regiones 
de alta densidad cerca del centro del cumulo, y por lo tanto la lumi- 
nosidad volumetrica podria estar infravalorada. Esta cuestion debe 
ser investigada en mayor detalle. 

13. Estimaciones observacionales tambien dan luminosidades menores a 
temperatura dada. La discrepancia es de un factor 2 para la esti- 
mation basada solo en la cosmologi'a ACDM, y es mucho peor para las 
estimaciones SCDM. Sin embargo, la comparacion entre las diferentes 
cosmologias no es trivial. 

6.3.3 Estrellas 

Hemos realizado 66 simulaciones numericas agrupadas en siete experimentos 
diferentes, con el fin de obtener una prediction teorica de la historia de la 
formation estelar cosmica. Comparamos nuestras predicciones con estima- 
ciones observacionales publicadas. Resumimos nuestros resultados: 

1. ACDM muestra el mejor acuerdo con las observaciones. El modelo 
CDM Estandar tiende a sobreestimar la densidad cosmica SFR, mien- 
tras que la BSI falla al forma suficientes estrellas a gran corrimiento 
hacia el rojo debido a su fuerza reducida a escalas pequenas. 

2. La fotoionizacion fija el comienzo de la formation estelar en el un- 
viverso, suprimiendo la formation estelar de objetos masivos a alto z. 
La retroalimentacion de las supernovas es muy efitiente para autoreg- 
ular el proceso de convertir gas en estrellas, previniendo asf que los 
halos mas masivos formen demasiadas estrellas en epocas tardi'as. 

3. Al contrario de las primeras estimaciones observacionales, pero de 
acuerdo con los datos mas recientes, el modelo ACDM predice que no 
existe decaimiento en el SFR cosmico para 2 < z < 5. La formation 
estelar parece ser un proceso gradual sin ninguna epoca caracterfstica. 
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4. Los volumenes comoviles en los experimentos son comparables con las 
observaciones. Por lo tanto se puede esperar un error apreciable en las 
medidas observacionales debidas a la variancia cosmica asociada a la 
estadi'stica de volumenes pequenos. Esto es lo mas importante cuando 
se quiere determinar el salto de pendiente de la densidad SFR para 
z ~ 1 en el momento actual. 

5. La historia de la formacion estelar es muy diferente en objetos que 
acaban en los centros de los cumulos de las galaxias de campo. Las 
galaxias de los centros de los cumulos tienen una poblacion estelar mu- 
cho mas vieja y su actividad se reduce drasticamente cuando pierden 
la reserva de gas frfo a partir de la que se forman las estrellas. 

6. En concordancia con las observaciones mas recientes, encontramos que 
los mecanismos que inhiben la formacion estelar en cumulos de galax- 
ias son tambien eficientes en entornos correspondientes a grupos mas 
pequenos (M ~ 1O 13 M ). 

Los efectos de selecion, extincion por polvo y conversion de luminosidades 
a tasas de formacion estelar son las principales fuentes de incertidumbre de 
las estimaciones observacionales. asf mismo el pequeno volumen estudiado 
por el HDF tambien podrfa introducir algiin tipo de sesgo estadi'stico en 
las observaciones a alto redshift. 

Las incertidumbres en los experimentos numericos provienen fundamen- 
talmente de la resolucion espacial y de la modelizacion de la formacion es- 
telar y de los procesos de realimentacion. El uso de codigos eulerianos con 
malla adaptativa sera muy util para conseguir un aumento significativo de 
la resolucion espacial en volumenes cosmologicos que sean estadi'sticamente 
significativos. 

El enriquecimiento qmmico y los procesos de fotoionizacion han sido 
tenidos en cuenta de forma fenomenologica en el codigo que hemos utilizado. 
Seria extremadamente interesante realizar un modelo autoconsistente de am- 
bos procesos simultaneamente, calculando la intensidad del fondo ionizante 
ultravioleta, no solo proveniente de quasares y AGN, sino tambien de las 
estrellas generadas en la simulacion en cada paso de tiempo, asf como la 
produccion de metales en las mismas y su tratamiento hidrodinamico como 
otra variable mas. 
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Appendix A 

Cluster Atlas 



Is this the real life ? Is this just fantasy ? 
Caught in a landslide, no escape from reality 
Open your eyes, look up to the skies and see! 

- Queen : Bohemian Rhapsody (1975) - 



'l)tOl)0!)Mlt this work, the sample of galaxy clusters described in Sec- 
tion 12,3.11 has been analysed in a statistical way. The present appendix is 
devoted to test whether the theoretical considerations that have been dis- 
cussed in the previous chapters can be applied to predict the structure of 
individual clusters. 



A.l Figure description 

We have been plotted two figures per cluster. The first one is divided into 
two columns, which in turn have been sub-divided into three panels. The 
second one contains a contour plot overlaid on a colour map. 

The left column of the first figure is devoted to the dark matter distri- 
bution, while the right column is used to display the radial properties of the 
ICM gas. In both cases, the numerical results are indicated by dots, while 
solid lines are used to represent the theoretical predictions. 

Dark matter: Dark matter density is plotted on the top panel, nor- 
malised to the mean cosmic value. In the middle panel, the cumulative 
mass is shown, and the lower panel is used to display its logarithmic slope 
a.M- The solid lines show the predictions based on the spherical collapse 
formalism, as discussed in Section [3.3.41 The values of the peak height and 
smoothing scale appropriate for each cluster can be found in Table 13.21 We 
recall that major mergers have not been fit by the spherical infall model. 

Gas: Top panel shows the gas density, also normalised to mean cosmic 
value. Solid line corresponds to the prediction 1)4. 17|) based on hydrostatic 
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equilibrium with a NFW-like potential, a polytropic equation of state and 
constant baryon fraction at large radii. Middle panel is used for the tem- 
perature, which has been compared with our prescription (|4.15|) . The lower 
panel shows the entropy profile of the cluster, computed as S = kTn e 2 ^ 3 
both for the experimental and theoretical profiles. As has been shown in 
Section 14.11 the entropy-conserving version of Gadget is able to produce 
isentropic cores in many of our clusters. For details on the computation 
of the theoretical predictions of the gas profiles, the reader is referred to 
Section E221 

Colour plot: This image represents the a map of the dark matter 
density, projected along the three axes. As can be appreciated in some fig- 
ures, the cluster might look more or less 'relaxed' (where, in this context, 
we understand 'relaxed' as 'spherically symmetric') depending on the pro- 
jection observed. If we assume that light traces mass, this colour plots can 
be thought of as a coarse approximation to an observational image of the 
cluster in an optical band. 

Contour plot: These contours correspond to the projected X-ray lu- 
minosity of the clusters. Although they approximately trace the observed 
shape of the dark mater distribution, it is evident that the gas has much 
less substructure than the dark matter. This is due to the gas stripping that 
galaxies experience as they orbit through the cluster. Most of the infalling 
haloes lose their gas reservoir during the first pericentre passage. 



A.2. INDIVIDUAL CLUSTERS 

A. 2 Individual clusters 
Cluster A 
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Figure A 
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Cluster A 
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Cluster B 
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Figure A. 2: Cluster B 
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Cluster C 
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Cluster D 
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Figure A. 4: Cluster D 
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Cluster E 
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Figure A. 5: Cluster E 
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Cluster F 
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Figure A. 6: Cluster F 
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Cluster G 
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Figure A. 7: Cluster G 
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Cluster H 
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Figure A. 8: Cluster H 
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A.2. INDIVIDUAL CLUSTERS 
Cluster I 
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Figure A. 9: Cluster I 
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Cluster Ji 
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A.2. INDIVIDUAL CLUSTERS 
Cluster Ki 
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A.2. INDIVIDUAL CLUSTERS 
Cluster L 
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Figure A. 12: Cluster L 
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Cluster M 
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Figure A. 13: Cluster M 



APPENDIX A. CLUSTER ATLAS 



A.2. INDIVIDUAL CLUSTERS 
Cluster K 2 
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A.2. INDIVIDUAL CLUSTERS 
Cluster J 2 
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